{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "\\newcommand{\\vt}[1]{\\boldsymbol{#1}}\n",
    "\\newcommand{\\uv}[1]{\\hat{\\boldsymbol{#1}}}\n",
    "\\newcommand{\\tv}[1]{\\tilde{\\boldsymbol{#1}}}\n",
    "\\newcommand{\\dvg}{\\operatorname{div}}\n",
    "\\newcommand{\\grd}{\\operatorname{grad}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\renewcommand{\\vt}[1]{\\boldsymbol{#1}}\n",
    "\\renewcommand{\\uv}[1]{\\hat{\\boldsymbol{#1}}}\n",
    "\\renewcommand{\\tv}[1]{\\tilde{\\boldsymbol{#1}}}\n",
    "\\renewcommand{\\dvg}{\\operatorname{div}}\n",
    "\\renewcommand{\\grd}{\\operatorname{grad}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integration of $R^n$ on triangles\n",
    "\n",
    "We are interested in the computation of\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_T R^n dy,\n",
    "\\end{equation}\n",
    "\n",
    "where $T$ is a triangle, $R=|x-y|$, and $n = -3,-1,0,1,2,3,..$. This package is based on [1] and papers building on this paper, notably [2,3].\n",
    "\n",
    "[[1]](http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1143304&tag=1) D. Wilton, S. Rao, A. Glisson, D. Schaubert, O. Al-Bundak, and C. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Transactions on Antennas and Propagation, vol. 32, no. 3, pp. 276–281, Mar. 1984.\n",
    "\n",
    "[2]R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D Green’s function or its gradient on a plane triangle,” IEEE Transactions on Antennas and Propagation, vol. 41, no. 10, pp. 1448–1455, Oct. 1993.\n",
    "\n",
    "[3] P. Yla-Oijala and M. Taskinen, “Calculation of CFIE impedance matrix elements with RWG and n times;RWG functions,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 8, pp. 1837–1846, Aug. 2003.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following Wilton[1], we start with writing the integrand as the surface divergence of a radially symmetric vector function. Assume\n",
    "\n",
    "\\begin{equation}\n",
    "R^n = \\dvg_{y} \\left( F(P;h) \\vt{P} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "Here, $\\vt{P}$ is the projection onto the plane supporting $T$ of $\\vt{y} - \\vt{x}$, and $h$ is the signed normal coordinate of $\\vt{x}-\\vt{y}$, i.e. $\\vt{x}-\\vt{y} = h\\uv{n} - \\vt{P}$ (with $\\vt{n}$ a normal on $T$ that will induce a counter-clockwise ordering of $\\partial T$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Ansatz above is equivalent to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\left(P^2+h^2\\right)^{n/2} = F'(P;h) P + F(P;h) 2\n",
    "\\end{equation}\n",
    "\n",
    "Multiplying with the integrating factor $P$ leads to\n",
    "\n",
    "\\begin{equation}\n",
    "P \\left(P^2+h^2\\right)^{n/2} = \\left( P^2 F(P;h) \\right)'\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or,\n",
    "\n",
    "\\begin{align}\n",
    "F(P;h) =& \\frac{1}{P^2} \\int Q \\left(Q^2+h^2\\right)^{n/2} dQ \\\\\n",
    "=& \\frac{1}{P^2} \\frac{1}{(n+2)}  \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}}\n",
    "\\end{align}\n",
    "\n",
    "and finally\n",
    "\n",
    "\\begin{equation}\n",
    "R^n = \\dvg \\left( \\frac{1}{P^2} \\frac{1}{(n+2)} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} \\vt{P} \\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation by Recursion\n",
    "\n",
    "Using this the integral we want to compute reduces to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\lim_{\\epsilon \\rightarrow 0} \\frac{1}{n+2} \\int_{\\partial(T \\setminus B_{\\epsilon})} \\frac{ \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}}}{P^2} \\uv{m} \\cdot \\vt{P} dy - \\frac{1}{n+2} \\lim_{\\epsilon \\rightarrow 0} \\int_{C_{\\epsilon}} \\frac{ \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}}}{P^2} \\uv{m} \\cdot \\vt{P} dy\n",
    "\\end{equation}\n",
    "\n",
    "The integration path of the contour in the second term is $\\partial (B_{\\epsilon}(\\xi) \\cap T)$. For $\\xi \\notin T$ the contour is void, for $\\xi \\in T$ it is a circle, for $\\xi$ in the interior of one of the sides comprising $\\partial T$ it is half a circle, and for $\\xi$ in the vertices of $T$ it is a circle spanning the angle formed at those vertices w.r.t the interior of $T$. If $\\xi \\in \\partial T$, the integral in the first term is meant in the Cauchy principle value sense.\n",
    "\n",
    "The second term is a contribution of a small circlular area of radius $\\epsilon$ centered on $\\xi$. This region needs to be excluded prior to applying the surfacic divergence theorem because the potential is not smooth there. It can easily be evaluated to\n",
    "\n",
    "\\begin{align}\n",
    "S_n =& - \\frac{1}{n+2} \\lim_{\\epsilon \\rightarrow 0} \\int_{C_{\\epsilon}} \\frac{ \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}}}{P^2} \\uv{m} \\cdot \\vt{P} dy \\\\\n",
    "=& -\\frac{\\alpha}{n+2} |h|^{n+2}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary of the triangle is the disjoint union of the three segments: $\\partial T = \\cup_{i=1}^3 \\partial_i T$. On $\\partial_i T$, $\\uv{m} \\cdot \\vt{P}$ takes on the constant value denoted $p_i$ (the sign of $p_i$ depends on the position of $\\vt{\\xi}$ with respect to the oriented line $\\partial_i T$). The segment $\\partial_i T$ is parametrised by\n",
    "\n",
    "\\begin{equation}\n",
    "\\partial_i T(s) = \\vt{y}_{0,i} + s \\uv{t}_i, \\quad s \\in [a_i,b_i]\n",
    "\\end{equation}\n",
    "\n",
    "where $\\vt{y}_{0,i}$ is the orthogonal projection of $\\vt{x}$ on $\\partial T_i$, and $\\uv{t}_i$ is the unit vector along $\\partial_i T$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integral reduces to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{n+2} \\sum_{i=1}^3 p_i \\int_{a_i}^{b_i} \\frac{\\left(s^2 + p_i^2 + h^2\\right)^{\\frac{n+2}{2}}}{s^2 + p_i^2} ds\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now focus on a single term:\n",
    "\n",
    "\\begin{equation}\n",
    "I_n = \\frac{p}{n+2} \\int_{a}^{b} \\frac{(s^2+p^2+h^2)^{\\frac{n+2}{2}}}{s^2+p^2} ds\n",
    "\\end{equation}\n",
    "\n",
    "where indices $i$ where omitted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $n = -3$ we get\n",
    "\n",
    "\\begin{align}\n",
    "I_{-3} =& -p \\int_{a}^{b} \\frac{(s^2+p^2+h^2)^{-\\frac{1}{2}}}{s^2+p^2} ds \\\\\n",
    "=& -\\frac{1}{h} \\tan^{-1} \\frac{hs}{p \\sqrt{s^2+p^2+h^2}}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note*: Here we used the dare I say non-trivial primitve:\n",
    "\n",
    "\\begin{equation}\n",
    "\\int \\frac{ab}{\\left( s^2 + a^2 \\right) \\sqrt{s^2 + a^2 + b^2}} = \\tan^{-1} \\frac{bs}{a \\sqrt{s^2 + a^2 + b^2}}\n",
    "\\end{equation}\n",
    "\n",
    "See e.g. [Wolfram Alpha](https://www.wolframalpha.com/input/?i=diff%28atan%28b*s%2F%28a*sqrt%28s%5E2%2Ba%5E2%2Bb%5E2%29%29%29,s%29).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For n = 0 the integral becomes\n",
    "\n",
    "\\begin{align}\n",
    "I_0 =& \\frac{1}{2} p \\int_a^b \\frac{s^2+p^2+h^2}{s^2+p^2} ds \\\\\n",
    "=& \\frac{1}{2} p \\left( s + \\frac{h^2}{p} \\tan^{-1} \\frac{s}{p} \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case $n > 0$ the integrals can be computed by recursion.\n",
    "\n",
    "\\begin{align}\n",
    "I_n =& \\frac{p}{n+2} \\int \\frac{(s^2+p^2+h^2)^{\\frac{n+2}{2}}}{s^2+p^2} ds \\\\\n",
    "=& \\frac{p}{n+2} \\int \\frac{(s^2+p^2+h^2)^{n/2}\\left(s^2+p^2+h^2\\right)}{s^2+p^2} ds \\\\\n",
    "=& \\frac{p}{n+2} \\left( \\int \\left( s^2+p^2+h^2 \\right)^{n/2} ds + h^2 \\frac{n}{p} I_{n-2}\\right) \\\\\n",
    "=& \\frac{p}{n+2} J_n(s) + h^2 \\frac{n}{n+2} I_{n-2} \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first term can be computed by recursion. Consider the sequence of integrals:\n",
    "\n",
    "\\begin{align}\n",
    "L_n =& \\int \\left( t^2 + 1 \\right)^{n/2} dt \\\\\n",
    "=& \\int \\cosh^{n+1} (\\alpha) d\\alpha \\quad \\mbox{where } t = \\sinh \\alpha \\\\\n",
    "=& \\cosh^n \\alpha \\sinh \\alpha - n \\int \\cosh^{n-1} \\alpha \\sinh^2 \\alpha \\\\\n",
    "=& \\left(1 + t^2\\right)^{\\frac{n}{2}} t - n \\int \\cosh^{n-1} \\alpha (\\cosh^2 \\alpha - 1) \\\\\n",
    "=& \\left(1 + t^2\\right)^{\\frac{n}{2}} t - n \\left( L_n - L_{n-2} \\right)\n",
    "\\end{align}\n",
    "\n",
    "From which\n",
    "\n",
    "\\begin{equation}\n",
    "L_n = \\frac{1}{n+1} \\left( \\left(1 + t^2\\right)^{\\frac{n}{2}} t + n J_{n-2} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "The recursion is bootstrapped from\n",
    "\n",
    "\\begin{equation}\n",
    "L_{-1} = \\log \\left( t + \\sqrt{1 + t^2} \\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And so we get for the sequence $J_n$:\n",
    "\n",
    "\\begin{align}\n",
    "J_{n}(s) =& \\int \\left( s^2 + p^2 + h^2 \\right)^{\\frac{n}{2}} ds \\\\\n",
    "=& \\left( p^2 + h^2 \\right)^{\\frac{n+1}{2}} \\int (t^2+1)^{n/2} dt, \\quad t = \\frac{s}{\\sqrt{p^2+h^2}} \\\\\n",
    "=& \\left( p^2 + h^2 \\right)^{\\frac{n+1}{2}} L_n(t) \\\\\n",
    "=& \\frac{1}{n+1} \\left( p^2 + h^2 \\right)^{\\frac{n+1}{2}} \\left( \\left(1 + t^2\\right)^{\\frac{n}{2}} t + n L_{n-2} \\right) \\\\\n",
    "=& \\frac{1}{n+1} \\left( s \\left( s^2 + p^2 + h^2 \\right)^{\\frac{n}{2}} + n (p^2+h^2) J_{n-2}(s) \\right)\n",
    "\\end{align}\n",
    "\n",
    "bootstrapped by\n",
    "\n",
    "\\begin{equation}\n",
    "J_{-1}(s) = \\log \\left( s + \\sqrt{s^2 + p^2 + h^2} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "This recursion was already mentioned in [3]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case where $p=0$ and $h=0$ we will use the limiting value\n",
    "\n",
    "\\begin{equation}\n",
    "\\left. J_{-1}(s) \\right|_a^b = \\lim_{\\epsilon \\rightarrow 0} \\log \\left( \\frac{b+\\sqrt{b^2+\\epsilon^2}}{a + \\sqrt{a^2 +\\epsilon^2}} \\right) = \\left\\{ \\begin{array}{cl}\n",
    "\\log \\left( \\frac{b}{a} \\right), & \\mbox{when } b > 0 \\\\\n",
    "\\log \\left( \\frac{a}{b} \\right), & \\mbox{when } b < 0 \\\\\n",
    "\\end{array} \\right.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*note*: The limit is only finite if $a$ and $b$ have the same sign. This is guaranteed because the singularity $\\vt{x}$ is not allowed to be on the boundary (even though it is allowed to be on the extension of these segments)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first recursion gives the integral for $n = -1$, i.e. one of the values computed in [1]:\n",
    "\n",
    "\\begin{align}\n",
    "I_{-1} =& p \\int_{a}^{b} \\frac{(s^2+p^2+h^2)^{1/2}}{s^2+p^2} ds \\\\\n",
    "=& p \\int_{a}^{b} \\left( \\frac{1}{\\sqrt{s^2+p^2+h^2}} + \\frac{h^2}{\\sqrt{s^2+p^2+h^2}(s^2+p^2)} \\right) ds \\\\\n",
    "=& p \\left. \\log (s + \\sqrt{s^2+p^2+h^2}) \\right|_a^b + h \\left. \\tan^{-1} \\frac{hs}{p \\sqrt{s^2+p^2+h^2}} \\right|_a^b\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combining the contour and singularity contributions\n",
    "\n",
    "For small $p$, the argument of $\\tan^{-1}$ becomes very large. Moreover, the contributions stemming from the integral over $\\partial T$ and those stemming from the region around $\\xi$ cancel in the final computation of $\\sum_i I_{n,i} + S_n$.\n",
    "\n",
    "In [1] it is suggested to combine both contributions using the addition property of the $\\tan^{-1}$. In this section this will be generalised for arbitrary powers of $R$. Not only will this result in the stable computation of the integrals, it will also eliminate branches from the code.\n",
    "\n",
    "Using the computations from previous section we have\n",
    "\n",
    "\\begin{align}\n",
    "\\int \\frac{h}{R^3} =& - \\sum_{i=1}^3 \\left. \\tan^{-1} \\frac{hs}{p_i \\sqrt{s^2+p_i^2+h^2}} \\right|_{a_i}^{b_i} + \\alpha \\sigma(h) \\\\\n",
    "=& - \\sigma(h) \\left( \\sum_{i=1}^3 \\left. \\tan^{-1} \\frac{|h|s}{p_i \\sqrt{s^2+p_i^2+h^2}} \\right|_{a_i}^{b_i} - \\alpha \\right) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interior angle spanned by $\\partial T$ as seen by $\\xi$ can be written as:\n",
    "\n",
    "\\begin{equation}\n",
    "a = \\sum_{i=1}^3 \\left. \\tan^{-1} \\frac{s}{p_i} \\right|_{a_i}^{b_i}\n",
    "\\end{equation}\n",
    "\n",
    "The [addition theorem](https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Arctangent_addition_formula) for $\\tan^{-1}$ applied to this situtation allows to combine corresponding terms from $\\sum_i I_{n,i}$ and $S_n$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\tan^{-1} \\frac{s}{p_i} - \\tan^{-1} \\frac{|h|s}{p_i \\sqrt{s^2+p_i^2+h^2}} = \\tan^{-1} \\frac{p_i s}{(p_i^2 + h^2) + |h| \\sqrt{s^2+p_i^2+h^2}}\n",
    "\\end{equation}\n",
    "\n",
    "Note that\n",
    "\n",
    "\\begin{equation}\n",
    "S_n = h^2 \\frac{n}{n+2} S_{n-2}\n",
    "\\end{equation}\n",
    "\n",
    "and comparing this with the second term in the recursion for $I_n$ ($n$ odd):\n",
    "\n",
    "\\begin{equation}\n",
    "I_{n,i} = ... + h^2 \\frac{n}{n+2} I_{n-2,i}\n",
    "\\end{equation}\n",
    "\n",
    "tells us that once the merging of the two contributions is done for $n=-3$, the recursion will generate the correct answers for all other powers of $R$.\n",
    "\n",
    "Whenever $p=0$, the $\\tan^{-1}$ will be set to zero. Not only will this result in more efficient code in that case, it will also void division by zero when in addition also $h=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $n = 0$ the situation is even better:\n",
    "\n",
    "\\begin{align}\n",
    "I_{0} =& \\frac{1}{2} \\sum_{i=1}^{3} \\left( p_i \\left. s \\right|_{a_i}^{b_i} + h^2 \\left. \\tan^{-1} \\frac{s}{p} \\right|_{a_i}^{b_i} \\right) -\\frac{1}{2} h^2 \\alpha \\\\\n",
    "=& \\frac{1}{2} \\sum_{i=1}^{3} p_i \\left. s \\right|_{a_i}^{b_i}\n",
    "\\end{align}\n",
    "\n",
    "Again because the second term in the recursion for $I_n$ and the recursion for the singularity contribution behave identically, the $\\tan^{-1}$ will cancel  in the computation of integrals of all powers of $R$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The vector case\n",
    "\n",
    "In this section we will compute integrals of the form\n",
    "\n",
    "\\begin{equation}\n",
    "\\int \\frac{\\vt{y} - \\vt{\\xi}}{R^n} dy, \\quad n=-3,-1,0,1,2,...\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we will look for scalar potentials such that their surfacic gradient is the above integrand:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\vt{y} - \\vt{\\xi}}{R^n} = \\grd_T F(P;h)\n",
    "\\end{equation}\n",
    "\n",
    "or\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{P} F'(P;h) \\vt{P} = \\left( P^2 + h^2 \\right)^{n/2} \\vt{P}\n",
    "\\end{equation}\n",
    "\n",
    "from which\n",
    "\n",
    "\\begin{align}\n",
    "F(p) =& \\int \\left(P^2 + h^2 \\right)^{n/2} dP \\\\\n",
    "=& \\frac{1}{2} \\int_{Q = P^2+h^2} Q^{n/2} dQ \\\\\n",
    "=& \\frac{1}{n+2} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}}.\n",
    "\\end{align}\n",
    "\n",
    "So we conclude\n",
    "\n",
    "\\begin{equation}\n",
    "R^{n} = \\grd_T \\left( \\frac{1}{n+2} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "The integral we are looking for thus can be written as\n",
    "\n",
    "\\begin{align}\n",
    "\\int_T \\frac{\\vt{y} - \\vt{\\xi}}{R^n} dy =& \\lim_{\\epsilon \\rightarrow 0} \\left( \\int_{\\partial T \\setminus B_{\\epsilon}} \\frac{\\uv{m}}{n+2} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} dy + \\int_{T \\cap C_{\\epsilon}} \\frac{\\uv{m}}{n+2} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} dy \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $n \\geq -1$ the second term clearly disappears when $\\epsilon \\rightarrow 0$. Also when $n=-3$ and $h \\neq 0$ the singular contribution disappers. It is only the case where $n=-3$ and $h = 0$ that needs a more careful analysis of the second term. If $\\xi \\notin \\partial T$, the integral vanishes because contributions on opposite sides of $\\vt{\\xi}$ cancel. When $\\xi \\in \\partial T$ this is not necessarily the case. However, we can safely assume this case will never occur because in that case the integral $\\int_T (\\vt{y}-\\vt{\\xi}) R^{-3}$ diverges anyway.\n",
    "\n",
    "In conclusion either the integral diverges or the second term is zero. We can focus or efforts on the first term."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like we did for the scalar case, the integral can be written as\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_T \\frac{\\vt{y} - \\vt{\\xi}}{R^n} dy = \\sum_{i=1}^{3} \\uv{m}_i \\left. K_n(s,p_i,h) \\right|_{a_i}^{b_i}\n",
    "\\end{equation}\n",
    "\n",
    "with\n",
    "\n",
    "\\begin{equation}\n",
    "K_n(s,p,h) = \\frac{1}{n+2} \\int_a^b \\left( s^2 + p^2 + h^2 \\right)^{\\frac{n+2}{2}} = \\frac{1}{n+2} J_{n+2}\\left(s,p,h\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can be easily expressed in terms of quantities we already encountered in our previous caluclations. In other words computation of the vector integrals comes at almost no computational overhead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extension to Triangle/Disk intersections\n",
    "\n",
    "In the Space-Time Galerkin discretisation of time domain retarded potential boundary integral equations, the computation of integrals over the following class of domains is required:\n",
    "\n",
    "\\begin{equation}\n",
    "D = T \\cap B(\\xi,Q)\n",
    "\\end{equation}\n",
    "\n",
    "where $T$ is a triangle and $B(\\xi,Q)$ is the disk with centre $\\xi$ (assumed to lie in the plane supporting $T$) and radius $Q$.\n",
    "\n",
    "When given a triangle by its three vertices and a disk by its center and radius, the construction of a boundary representation of such a domain is nothing short of a nightmare, and will not be discussed in this document.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start from the alternative surfacic potential representation for the kernel:\n",
    "\n",
    "\\begin{equation}\n",
    "R^n = \\dvg \\left( \\frac{1}{(n+2) P^2} \\left( \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} - \\left(h^2 \\right)^{\\frac{n+2}{2}} \\right) \\vt{P} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "Note that the term we added has vanishing divergence and thus will not affect the value of the integral. The reason for this alternative potential is that it is regular at $P=0$ and so we do not have to keep track of the singularity term separately.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary $\\partial D$ is made up from straight segments and arcs. Integration of the potential along the straight segments is dealt with as explained in the previous sections. For an arc $C$ the contribution to the integral looks like:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "A_n =& \\frac{1}{n+2} \\int_{C} \\frac{\\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} - \\left( h^2 \\right)^{\\frac{n+2}{2}}}{P^2} \\uv{m} \\cdot \\vt{P} dy \\\\\n",
    "=& \\frac{\\alpha}{n+2} \\left( \\left( p^2 + h^2 \\right)^{\\frac{n+2}{2}} - \\left(h^2\\right)^{\\frac{n+2}{2}} \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\alpha$ is the signed angle spanned by the arc (positive for arcs traversed in counterclockwise sense, negative otherwise). From this expression we get the following recursion:\n",
    "\n",
    "\\begin{equation}\n",
    "A_n = \\frac{\\alpha}{n+2} \\left[ \\left( \\frac{n}{\\alpha} A_{n-2} + \\left( h^2 \\right)^{\\frac{n}{2}} \\right) \\left( p^2 + h^2 \\right) - \\left( h^2 \\right)^{\\frac{n+2}{2}} \\right]\n",
    "\\end{equation}\n",
    "\n",
    "Bootstrapped by\n",
    "\n",
    "\\begin{align}\n",
    "A_{-3} =& -\\alpha \\left( \\left( p^2 + h^2 \\right)^{-\\frac{1}{2}} - \\left( h^2 \\right)^{-\\frac{1}{2}} \\right) \\\\\n",
    "A_{0} =& \\frac{\\alpha}{2} p^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the vector case the contribution from a single arc reads\n",
    "\n",
    "\\begin{align}\n",
    "B_n =& \\int_C \\uv{m} \\left( \\frac{1}{n+2} \\left( P^2 + h^2 \\right)^{\\frac{n+2}{2}} \\right) dy \\\\\n",
    "=& \\frac{1}{n+2} \\left( p^2 + h^2 \\right)^{\\frac{n+2}{2}} \\int_C \\uv{m} dy\n",
    "\\end{align}\n",
    "\n",
    "*Note*: In the vector case we need to be consistent in our choice for the potential $\\Phi$ since the quantity we integrate over is $\\uv{m} \\Phi$ which does not necessarily cancel at the inter-sector boundary. In contrast, for the scalar case, the quantity $\\uv{m} \\cdot (F(P) \\vt{P})$ vanishes along the inter-sector boundary and does $F(P)$ can be chosen freely in each sector. So if we would leave out the second term in the surfacic potential, we would have to do so also in the computation of the segment contributions, forcing us to consider what happens in the neighborhood of $\\xi$. For certain configurations of $D$ this point might not even be in the integration domain. The resulting algorithm becomes inelegant and prone to numerical and logical errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note*: for contributions of arcs and circles we can assume that $p \\neq 0$. This ensures for example that no divisions by zero occur in the computation of $h A_{-3}$\n",
    "\n",
    "*Note* In the case of a full circle $\\alpha = \\pm 2\\pi$, $B_n$ vanishes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using WiltonInts84\n",
    "using FixedSizeArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p1 = Vec(0.0, 0.0, 0.0)\n",
    "p2 = Vec(1.0, 0.0, 0.0)\n",
    "p3 = Vec(0.0, 1.0, 0.0)\n",
    "h = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.5:0.01:1.5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = -0.5 : 0.01 : 1.5\n",
    "y = -0.5 : 0.01 : 1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0,0.15)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t, T = 0.15, 1\n",
    "r, R = max(0,t-T), max(0,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M = zeros(typeof(t), length(x), length(y))\n",
    "for i in eachindex(x)\n",
    "    for j in eachindex(y)\n",
    "        c = Vec(x[i],y[j],h)\n",
    "        try\n",
    "            I, K = wiltonints(p1,p2,p3,c,r,R,Val{0})\n",
    "            M[i,j] = I[2]\n",
    "        catch\n",
    "            #@show i j\n",
    "            #@show c\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnsAAAITCAYAAAB/k9qgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzsvXuQHXWd/v/05dwvEyAEhgj4hRBLDSESyEZuWRSCRHapxUSUFXErhipdSwwqIfpFq1BUFjCuUKULtWQ2wCosi1W73+j6K3QDRMhG1KBELrpCxGQSjZnMnPs5ffn90f0559N9+lx6Micz03leVR/Ome4+N+Z0+pn37VFs27ZBCCGEEEIiiTrdb4AQQgghhAwOij1CCCGEkAhDsUcIIYQQEmEo9gghhBBCIgzFHiGEEEJIhKHYI4QQQgiJMBR7hBBCCCERhmKPEEIIISTCUOwRQgghhEQYij1CCCGEkAgzI8RevV7Hhg0bMH/+fKTTaSxfvhxPPvnkdL8tQgghhJBZz4wQezfccAO+8Y1v4Prrr8c3v/lN6LqOVatW4dlnn53ut0YIIYQQMqtRbNu2p/MN7Ny5E8uXL8c999yD9evXAwBqtRoWLVqEk046Cdu3b5/Ot0cIIYQQMquZ9sje448/Dl3XsW7duua2RCKBtWvX4rnnnsPevXun8d0RQgghhMxupl3s7dq1CwsXLkQ2m/VsX7ZsWXM/IYQQQgiZHPp0v4HR0VEMDw+3bR8eHoZt29i3b1/g4w4ePIgf/vCHePOb34xUKjXot0kIIYQQAJVKBa+//jquuOIKzJ0796i//u9//3scPHhwYM8/d+5cnHbaaQN7/ulg2sVepVJBIpFo255MJpv7g/jhD3+ID33oQwN9b4QQQggJ5uGHH8bf/u3fHtXX/P3vf48Fp5+OxgBfI51O46WXXoqU4Jt2sZdKpVCr1dq2V6vV5v4g3vzmN7v3rgHwPID3DOT9kenmv8DfbVTh7zba8PcbXf4DwH7pOnz0OHjwIBpwrvyDiCkeBPBEuYyDBw9S7E0lw8PDgana0dFRAMApp5wS+LiWCHwewJ8BbJP2LgJw9hS+SzJ9JAG0p/lJFODvNtrw9xsNfgXgRd+2cQCdgzFHg7kAgtUBCWLaxd6SJUuwbds2FItFT5PGjh07oCgKlixZ0uMZ3gNH6H1wgO+SEEIIORY5G+3Bk80Afj8N76WFjsEImGkXRQNi2rtxV69eDcMwcP/99ze31et1jIyMYPny5Zg/f/40vjtCCCGEkNnNtIvYZcuWYc2aNdi4cSMOHDiABQsWYGRkBHv27MHmzZun++0RQgghZIahA4gN6HmjyIz4XA899BBuu+02PPzwwxgbG8PixYuxdetWXHjhhX0+w6KBvj8ynfB3G134u402/P1GlwWY7jQuCceMEHvxeBx33nkn7rzzzkk+A5sxogt/t9GFv9tow99vdFkA4MfT+g40DEbAaAN4zpnAjBB7hBBCCCH9wjRuOKa9QYMQQgghhAyOqIpYQgghhEQUpnHDwcgeIYQQQkiEYWSPEEIIIbMK1uyFg5E9QgghhJAIE1URSwghhJCIwpq9cDCyRwghhBASYRjZI4QQQsisgjV74WBkjxBCCCEkwkRVxBJCCCEkougYjICJqihiZI8QQgghJMJEVcQSQgghJKKwZi8cjOwRQgghhESYqIpYQgghhEQURvbCEdXPRQghhJCIwqHK4WAalxBCCCEkwjCyRwghhJBZBdO44WBkjxBCCCEkwkRVxBJCCCEkorBmLxyM7BFCCCGERBhG9gghhBAyq2DNXjgY2SOEEEIIiTBRFbGEEEIIiSis2QsHI3uEEEIIIRGGkT1CCCGEzCpYsxcORvYIIYQQQiJMVEUsIYQQQiKKjsEImKiKIkb2CCGEEEIiTFRFLCGEEEIiCmv2whHVz0UIIYSQiMLRK+FgGpcQQgghJMIwskcIIYSQWQXTuOFgZI8QQgghJMJEVcQSQgghJKIwshcORvYIIYQQQiJMVEUsIYQQQiIKu3HDwcgeIYQQQkif1Ot1bNiwAfPnz0c6ncby5cvx5JNP9nzcj3/8Y6xduxZvectbkMlkcOaZZ2LdunXYv39/27G2bePb3/423vGOdyCXy+Hkk0/GqlWr8Nxzz03qPVPsEUIIIWRWoWtATJ/6pfcR2rvhhhvwjW98A9dffz2++c1vQtd1rFq1Cs8++2zXx23YsAFPPfUUrrnmGtx777344Ac/iMceewznnnsu/vjHP3qO/cxnPoOPf/zjOOecc7Bp0yZ85jOfwauvvooVK1bg+eefD///K/QjCCGEEEKOQXbu3IlHH30U99xzD9avXw8AuP7667Fo0SLccsst2L59e8fHbtq0CRdddJFn2xVXXIEVK1bgvvvuw+233w4AME0T3/72t/H+978fIyMjzWNXr16NM844A4888gjOO++8UO+bkT1CCCGEzCo0DdD1qV9aj8je448/Dl3XsW7duua2RCKBtWvX4rnnnsPevXs7PtYv9ADg4osvxvHHH4+XXnqpua3RaKBSqWDevHmeY0888USoqop0Ot3n/6UWFHuEEEIIIX2wa9cuLFy4ENls1rN92bJlzf1hKJVKKBaLmDt3bnNbMpnEX/zFX2BkZAT/+q//ijfeeAO//OUv8ZGPfAQnnHCCR2j2C9O4hBBCCJlV6CoQG0DrbC9RNDo6iuHh4bbtw8PDsG0b+/btC/V6mzZtQqPRwAc+8AHP9kceeQTvf//78aEPfai57cwzz8T27dvx5je/OdRrAIzsEUIIIYT0RaVSQSKRaNueTCab+/vl6aefxu23345rr70WK1as8OzLZrN4+9vfjk984hP43ve+h29961swDANXX301Dh06FPp9M7JHCCGEkFmF3mfnbOjnVbrvT6VSqNVqbdur1Wpzfz+8/PLLuOaaa7B48WI88MADnn2maeKyyy7DpZdein/8x39sbn/3u9+Nt7/97bjrrrvw1a9+ta/XEVDsEUIIIWRWIUavHAnfqTlLZtzu/pjh4eHAVO3o6CgA4JRTTun5um+88QZWrlyJ4447Dlu3bkUmk/Hsf/rpp/Hiiy9i06ZNnu0LFizAW9/6VvzkJz/p+Rp+KPYIIYQQcszxwYSzZH5uAEvHOz9myZIl2LZtG4rFoqdJY8eOHVAUBUuWLOn6mocOHcLKlSthGAa2bduGk046qe2YAwcOQFEUmKbZtq/RaMAwjO4fLADW7BFCCCFkdqHC8Tab6tVDFa1evRqGYeD+++9vbqvX6xgZGcHy5csxf/58AMD+/fvxyiuveARbuVzGlVdeidHRUXz/+9/HGWecEfgaCxcuhG3b+O53v+vZ/vOf/xyvvPIKzj333J7/e/wwskcIIYQQ0gfLli3DmjVrsHHjRhw4cAALFizAyMgI9uzZg82bNzePu/XWW7Flyxa8/vrrOO200wAA1113HX76059i7dq12L17N3bv3t08PpvN4uqrrwYAnHvuubj88svxL//yLxgfH8fKlSuxb98+3HfffchkMrjppptCv2+KPUIIIYTMLjQMRsFYvQ956KGHcNttt+Hhhx/G2NgYFi9ejK1bt+LCCy9sHqMoClTVGyZ84YUXoCgKHnzwQTz44IOefaeffnpT7AHAf/zHf+Duu+/Gd7/7Xfzwhz9EPB7HJZdcgttvvx1nnXVW6I+l2LbdoxxxZvLzn/8cS5cuBXAjgPaZN4QQQggZBKMA7sfPfvazSaUUjwRx7f/ZXODc+ACevw4sPYhp+WyDhJE9QgghhMwudExbZG82wgYNQgghhJAIw8geIYQQQmYXg6rZa592EgkY2SOEEEIIiTCM7BFCCCFkdiHm7A3ieSNIRD8WIYQQQggBGNkjhBBCyGxjUDV7g4gWzgAY2SOEEEIIiTCM7BFCCCFkdjGoOXsRVUUR/ViEEEIIiSxs0AhFRD8WIYQQQggBGNkjhBBCyGyDDRqhYGSPEEIIISTCMLJHCCGEkNkFGzRCwcgeIYQQQkiEiaiGJYQQQkhkYTduKCL6sQghhBBCCMDIHiGEEEJmG+zGDQUje4QQQgghEYaRPUIIIYTMLhjZCwUje4QQQgghEYaRPUIIIYTMLjhnLxSM7BFCCCGERJiIalhCCCGERBbO2QsFxR4hhBBCZhds0AhFRDUsIYQQQggBGNkjhBBCyGyDkb1QMLJHCCGEEBJhGNkjhBBCyOxCw2CicIzsEUIIIYSQ2QYje4QQQgiZXbBmLxSM7BFCCCGERBhG9gghhBAyu2BkLxSM7BFCCCGERBhG9gghhBAyu2A3bigo9sgMQ3EXpNugfd22EUIGi+2uTtu67SOEHG0o9sgMQ4FTXaCgXcip8FYeBB1DCBk8NtoFnOWuTvso9sgUwpq9UFDskWnCH70Tok0IOr/oA5yz0C8ExXFBz00IOTKCInRAsLCTt4n7AnGuytso/gg5WlDskWnAL9j84k5FsOgTf8oFPUYQtI0QEg4bwRE68bPpriCB5xeCFrznurydgo9MEkb2QkGxR44y4h99De3CrtfS3RW0TyCLRELI5JFFG+AVdCYAA17BF3YpoNgjk4ZiLxQUe+QoIf9lL9Kx8tKk5T9GbNd9x8lpXX8kUH5NQkh4RPTOn6aVhZ4hHWdK+23fz/Jju52TbOQgZBBQ7JGjiCz2ZKEmRJwu3ZeFnn+/7nsOVdqm+F6LYo+QyREUuZNFnn+J7Zp7rOb+rMAbzTMQLOgo9EgIOHolFBR7ZMDIaVu/ePNH7XQAMd9+vcOKoT3KJ1K84nX9KV5CSHfkOj1ZxImoXAPBQk9s1+GN9Km+55Dx1/DJzR8UfYRMJRR7ZIAEpW3lKJ1f6OkA4mgJOVnYydti7nH+5+pVz0cI6Yy/+UIWe+LnBoC6tL3hLg3tUT9xX/w7YPpez9/Y4W8E8YtDQiRYsxcKij0yIOSIXlAErlO0Tog9WeTJy79del5FAxRX4DVv3ffBbC4hwdjuf2wAthvBs93VJvbqcM47OZInlhCD4pxsoPUHlxB9An9q1z+WxS8MCSFHAsUeGQD+2jx/PZ4/YievOFqCT94m75Ofy30tTQE0FVAV30L7WD5CiENbk6wCmBpgKoClAZYOb3NFHI7ga0irDud8FEJPiDwFrVRupxPPhjedS0ifMLIXCoo9MgD8TRjiZxGZk1O2srjzizr554S7kghM0crlff4scdDYPkKOdURZnglH3BkADMXRauKKZwOw5Vq6OlqCTyzdvfU3XolmDHHCdbJYk/cxdUvIIKDYI1OI+Os9qJM2Jt36hV0i4FY+xl1KDFBVZ2mKV8jJh8rBxBjQVsJHsUeOZYJ6MOQ+C6HhxM/NfgsFsHXAlks0dAA1dJ+PKQQh3Mc1fG9ING/Y7vEWugtEQjCt3bj1eh233XYbHn74YYyNjWHx4sX48pe/jMsuu6zr43784x/jkUcewfbt2/GHP/wBJ598Mt71rnfhS1/6Ek4++eSOjxsfH8dZZ52FgwcP4vHHH8c111wT9lNR7JGpxG955q/PE2pMROn8UTsRuROiT67T0520rKa2l/bp0kPk6J5f/FHsEeLtiRCZWL/Qq8PRcHVpO9zUrinObROdm6LkmgkRrZNTuUECTog81Xccu3PJzOKGG27AE088gfXr12PBggUYGRnBqlWrsG3bNlxwwQUdH7dhwwaMjY1hzZo1OOuss/C73/0O9957L7Zu3Ypdu3Zh3rx5gY+77bbbUK1WoSiTr0Oi2CNTRJDIEyorqCZPjuQJgZeUlhB4mrM0zSvsEvBmesW2oMbeoJSueMuEHCvI2Vh/c22Q0Iu7tzWldf6INK+lAJZci+cfbC62++3S/OlaOcyoSrcCy/dchLhMU83ezp078eijj+Kee+7B+vXrAQDXX389Fi1ahFtuuQXbt2/v+NhNmzbhoosu8my74oorsGLFCtx33324/fbb2x7z4osv4tvf/ja++MUv4gtf+EL4z+NCsUemAPEPvN/pwh/VkxWaEHspeIWee19xO2p1BYj7onnyoX7dGKQrg4w32KxBjjXk1K2I6InoXZDQ8/dO1dDK2DbcGj9Lg3PyybV6/heUBZsIKerSPjki6Bd78nMxwkemn8cffxy6rmPdunXNbYlEAmvXrsXnP/957N27F/Pnzw98rF/oAcDFF1+M448/Hi+99FLgY2666Sa8733vw0UXXQTbnvz3n2KPHCFBQ5ODxqvIikzcT8IRe0K5CaGX8NbiJdsPadvmj/TJL9etaYOCj0QZ+dogp25rvlWX7svnkf80bprUKE50z1YAW5xc/r+ggqJ4/vZf0a0bJPbkTl166RIf0xTZ27VrFxYuXIhsNuvZvmzZsub+TmIviFKphGKxiLlz57bt+7d/+zfs2LEDL7/8Mn73u9/1/ZxBUOyRKcBfqyencv3jVGTFlnaXq8qUuJuyRSuaIPSgPwDoX2JfUJZYvmgF2eoSEmWE5hK1eXIEryrdVtG50Uk2tZH/tpPnKzdP3CBB5xd6shCUI3ziVgg9CjwysxgdHcXw8HDb9uHhYdi2jX379oV6vk2bNqHRaOADH/iAZ3u1WsVnP/tZ3HzzzTj11FMp9sh04xd6/nRu0EgVod7SaKm4uDMUWVO8qVr5sJRvu1/wSU/lifb5+jza0rmERBmhr2SRJ5Ys8uLufbnOtVM0XC7Ns+F06jZDIn4RJ4s+q8M+seQonnguQgKYpshepVJBIpFo255MJpv7++Xpp5/G7bffjmuvvRYrVqzw7PvqV78KwzCwcePGvp+vGwMTe0899RQuvfTStu2KouC5555rhjwB4OWXX8anPvUp/OQnP0E8Hsd73/tefP3rXw8Ma5KZiBzR6zQ42d9tK+VglZjTZaurrcOEmEtLK9Vh+cVfW2rXAuIWFN2GottQdec+NNsx2IC4JSQa2LbzhbYtBZapwjJU2DUVqKstgVdVgksggoRex5I8xdVudmsQsyfC5xd48goSe/J2VfqZRbbExzSNXkmlUqjVam3bq9Vqc38/vPzyy7jmmmuwePFiPPDAA559r7/+Ou6++25861vfQjqd7u9992Dgkb1PfepTOO+88zzbFixY0Ly/d+9eXHzxxTjuuOPwta99DYVCAXfddRdefPFF7Ny5E7rO4OPMRo7sBdXryRG9oPEqcXeciuot40ujXeh1u++P7skvkbCgxE2oMQNazIQWM6FqFlTNgqLaUBRGD0j0sG0FpqnBqOswGjrsmg6rpgNVFSjDEXwVtKYg+TvX5fv+BtugrGzD7dS19Q4H+aN84kWConua9BwCDlwmU8t3djlLZrza/THDw8OBqdrR0VEAwCmnnNLzdd944w2sXLkSxx13HLZu3YpMJuPZ/4UvfAFvetObcMkll2DPnj2e5//Tn/6EPXv24LTTTgs1imXgSuqiiy7qOgDwjjvuQKVS8RQ1nn/++bj88ssxMjKCj370o4N+i+SI6Geunl/ouUpMcVO3uiT0MmgJuX7upwGkbTeyZwMpG0rShpq0ocRtKHELatyAGmtAizWg6wZ03YCmmVBVE6pqQ1UsCj4SCWxXkVm2CstWYRoa6o04Go0YzHocRiMOq6LDjiuwyyoQV9wmJiW4a71Tg62s42y4zRpwonumjVa3rZjIbPluDfcFLN8Lysep0ovJfoc8VwmmJI37wfOcJfPzPwBLN3V+zJIlS7Bt2zYUi0VPk8aOHTugKAqWLFnS9TUPHTqElStXwjAMbNu2DSeddFLbMW+88QZ++9vf4owzzvBsVxQFH/vYx6AoCsbGxpDP53t/SJejEjYrFotIpVLQtPb46BNPPIGrrrrK073y7ne/GwsXLsRjjz1GsTdjkSN6cmTPH82T6/V8RXaqW6Mn0rZCyPVa6YBjUxaQsqCmTOhJw1musNO1OnStjpjaQExtQFcb0FQLmmJCVSyocBYhsxkbCiwosKHChAZT0WBoOuqIo64lUI8lUDcTaMTjaMR0mLEY7LjmRNWDmjBkodctSCfvM4Dm8GVbiD4drQHMwj6tn+VP38qduRR8ZHpYvXo17r77btx///24+eabATiOGiMjI1i+fHlTy+zfvx/j4+NYsGBBU/uUy2VceeWVGB0dxbZt29rEnOCOO+7AwYMHPdtefPFF3HbbbdiwYQPe+c53tkUDezFwsfd3f/d3KBQK0DQNF198Me666y4sXboUALBv3z788Y9/bEvzAk4b8w9+8INBvz0yafzdt/40rn/YnS+FqyRapT1yc24GQFa69d/Pwonk+Y9NW0DahJqqI5asIx6vIaE5K67UkUANcdQRQwNx1N3LoSUui1B58SCzHCeOpsGCBkPRYUBHQ4uhpiVQQwJVJJ1VS8KOJWG5p6itq4Bmu0JPaf87zh/F6zRBRczvs1UpAyvEnoHOhtV+5w2/4BPIxtYUfMc809SgsWzZMqxZswYbN27EgQMHmg4ae/bswebNm5vH3XrrrdiyZQtef/11nHbaaQCA6667Dj/96U+xdu1a7N69G7t3724en81mcfXVVwNAoAvH0NAQbNvG+eefj7/+678O/bEGJvbi8ThWr16NVatWYe7cufj1r3+Nu+++G5dccgmeffZZnHPOOc0cdKc25kOHDqHRaCAWiw3qbZIjIkjwyaIvyAdXd4Ylq/COVvFH63LwCjyxcmiKPCVjQcsaUDMG9EQdeqKGRLyOhF51hJ5ag3OZq3nEXgwN6DCgNYWeCRWWe43jBYTMLmwoblRPheF+sw33m15HvHkGNL//WgN63EANBhpqHA0tAVvTYalucZ4/VWv6lvDRlX9uoHW6y6P1PE1b/n8fhOVaJxcOQmYmDz30UJs37tatW3HhhRc2j1EUBarqHRD+wgsvQFEUPPjgg3jwwQc9+04//fSm2OvEjLRLe+c734l3vvOdzZ+vuuoqvO9978PixYuxceNGfP/732+2KPdqY6bYm4nINTRBgi+oQcOtAO8k9oTQk4WdJO6QA5Bv3VdyFrRMHfFMBUm9iqReQVKtOkupIoEakqh2FHtewWc1L5mEzCbs5p8rmvQNjzfFnojoiT0xrQE94ZQ4VLQUbN2GqSZgK04KOHBiioHuQq+Olthr6jd/45a/tlccGJSqZfct6cE0deMCTjDrzjvvxJ133tnxmM2bN3sifQDw2muvTfptrVixAqZpTvrxR7XV9cwzz8TVV1+N733ve7Btu9miPBVtzGQ66DZM2Z/Kde+rmmuBBu/AZDmyJ4u9puizXbFnQ82Z0HIGYpk6EqkykskyUmoZKbWClFJByr28pVBpE3udBF9L7DGyR2YHLatbzf0Wu/V50qoh0f59V5xudKcpyQZUoG4DsFWYlgLbVABT9Qo6f1RPLL/Qi6HVe2E6z+ktCNTdB/lTuUrAIoRMFUd9rsmpp56Ker2OUqnUTN+KdK7M6Ogojj/++D6iev8FRzXILAJw9lS8XdIT/+iVLkJPiL24EjxixR/ZE+ncvLtyNpC3oOfqiOcqSCUrSOklpNQy0koFKZSRQsWz0qggjrpP8NUQcxNd4gIoxB4hswVHqikw3W9xAzHUkJCEXrwZ0ZOj2M3vuZiDrACK5XbwWgBMHbYs9oKieMJb1++fG5OOMxXAUp1O3baJ5vIfh0zdzmx+BeBF37Ye80mOBtNUszdbOepi73//93+RTCaRzWaRzWZx4okn4vnnn287bufOnT1bmB3eA6C95o9MPbEYkMkAmYyCUklBqQQ0GkFiT/6XXzRnuP5LmtZug9ZN8OVsV+zZUPImlJyJWLaKdLaIdLyEDJyVRgXpptgrN++nUfHU7TmCr+YmuRqwKjbqJQCGjVS2gXSmwQHLZMbjr9PzR/OchoyEK/TkUgWnscGZh6zA1hRAA2xLhWFrsCwVpqHANrSWaDOU9kheUERPrAZaPRm2qFnqZsfBiN7M5my0B09GAdw/De+FTJaBib2DBw+2OWC88MIL+M///E+8973vbW573/vehy1btmDv3r3NluUf/ehHePXVV/HpT396UG+PTII5c4DlyzNYvvxk7Nx5EDt2NHDggD+iF2SV5rpkqO4sL7k5V3bCCErjuilcJW8inqsinqsikygioxaRQRFZlJBBEWmUkXFFniP0yk0B2OpFbNXwia7cPa9n8asdJ6BwSMc7lh/AO/7iIDSddXtkpuIbrwK1WZtXl77pTrlC3RO1lmtSbajucrYamg4jrsOyNNgNDVbDcockKy2R14BjsZaAI/Ti8Ao+v5bzaLcgW0UKO3IEMLIXioGJvWuvvRapVAoXXHAB5s2bh927d+OBBx5ANpvFV7/61eZxn/vc5/D444/jL//yL3HTTTehUCjg7rvvxjnnnIOPfOQjg3p7ZBIcd5yGiy8+BevWvRep1P+H3/zmDRw44Bd6/vtuGE9VXUs0eMetdGrQaKZvndStkjOQyFWRyRWQUYvIKkVkUUQWBVfwiQhfyRV9TrTPqduriKET7mo1a7zyWgy//M809r2WxqmpUcw9/0+IsUmDzFCctK3q1unpblQv5kbzks2u21b61nQlHaRSBUcsWu4WE6oj9jQdBmIwGzEYjRjQUFuuGMJXVwi9GoKjejqkUjwFaA4r9ws9pm0JOZoMTOz9zd/8DR555BFs2rQJExMTOPHEE7F69Wp84Qtf8AwSfNOb3oSnnnoKN998MzZu3Ih4PI6rrroKd999N7twZxjFooUXXzyI//f/tuNXvzqIQsH/D3mQe4b7r7+qtos8Wex50rl2s17PacaoI5GtIO1G9HJqATkUXLHnrAxKzVsh+kRkLwm3acN2RZ9dQ9w0EDMaWJhT8JdvV3BobgKLThjFUKEEnZE9MtNwA2G2osBSVdiqioamoaFqqKtxVFFHXPGOFRLizm6KKtsViy2hKO43lBjiiKGhxmDEdDQSOuyEDrumADEl2De329i8ZuCOgo4MCJFUGsTzRpCBib1PfOIT+MQnPtHXsW9961s5QHkWMDZm45lnxvGb3/wSBw4o+POfdXjr9cStPG7FFXua0krfygIvqG4vrQBZG8gCWs5AIldBOlNERnMiejm0iz3/yqDk1uy16viSqCJlV5G0aojXTcSqJs6eW8HwijHUyypOnFfB0EQZKq3TyEzDPbWc5lYFtq6gEVdYAP+uAAAgAElEQVRRj2uoKXHElHrXhiMh+pwaP82zGm5zR8MVew09hno8BiuuwIzp7ZE7v61ap3nIzNQSMmM46g0aZPZSqQCvvWbhtddEwY5/1p58BfD5L8ka0G+T2zZU2YaSsaDkLMQydaRSFaQTJbdGr7X8os8jAO1is34vaVWRsqpImjXEzRoSRh2xig29bCGrVHHySQW3mBzAoaP2v5OQ3sinlxQ0t2OAnlSgWyr0mAlNM6FrThWfpphNr2cR1bOhuOJObw5cFquOeKu9Q6mjpjtDl42YDitmw5ZFXieh103w9fyAYZAnPpNjGvF9HMTzRpCIfixydAiqw5EbNty5DkDrYhXknhbQoKFlDWjuHL2UVnLTtKUAcVdoir58QGo3Y5eQbNSRqDUQqxnQaia0qgWtAihltGqRTPD6QWYe/lNKlEIkALVuI1a3oCQaUBIWtLgJRbOhaK0onlgGdCRQRwO1jkIvhgZiijt0OWbAilnt0Ty/h64cwQsqyeup5cKcdGLaM09UAjZohIRij0wS/7iEILs0cUVAS+zJ0b22yJ7dbNJQMw3EsxUkk2Wk1bJbi1dsRveE6HNE3kRT7OVsSfBZzorXLMRLJrQSgDIAcVuBU2heh2vgTsgMQszB8ze3JwElCWgJQKvbUNNORE9XDSeip9mwbbf3VlFbdXk+kZdAHTWP50bdFXtONZ8Rs6DEbNgx2xV5Su+oniz0Jl371EnMeXzYCCEhoNgjk8SEo5J0OGoNaG/S0JxaPRXeaJ68UvA1alhA2oKeqCOpVR1nDMWZl5dxx6uIJgw5lZsXos8uImsVkDVLSNbqiNcs6EULShFAEd3FHrNDZCbQIXXbPIfEGBTX2UKx3alGio0YTCS0OkzVSdcmpEhe3WnB8A0UN9Fq1bCgwXLM1xQTqmIFj8Hr9HfepGr0/A/we7VZYESPBMLIXigo9sgkMeH8Ayw8LgFvKEIDFK01bsWfvvV35aYApGxX7JmIJepI6hXHAk2anydEXiuqJwk9N7KXMwvI1CqIFV2hVwCUCThiTwi9ILHH6wmZCfgjejq854vPn1aBc5opqo2YZsCO2zAVDYbSctaQ07W61LHrSMKWu4ZT82c5dX+qBah2sPV1J6ez0PhPuG5CjycnIZOFYo9MEn8Xrn9pgKK2Bin7jTWC0rlJG2rKhJaqIx6vIalWm1E90VkrR/Waos8uIGc7t9l6CalqDfGyAW0CUCdcoSeLPVnwCbHXAK8nZPrxR/TEyBMRyZP/KBGuZ+7xig7obp2dETPQ0BuoK3LyttFsz2j5QnuFnqo4kT1Fcb1zVdsJHSpKb7MLO2B13Wl1+NkK2McTk/jg6JVQUOyRSSK314pJquJff/dKIMSeXOcTYK4hBJ+StKGnDMSSNSR0V+yhiqQr9NLusGS/2HPq85yIXqpaQ2LCgFaQRJ5YBThCT47wVeGt2eM1hUwXQkD5z484Wv60lu9YqRZW0QEtZkPVLRgw0FAbiKlyJM9oirzWMiUvDtHO4XprKBYUxY3uNV80gKDMq6zX2g4MEnaM4hEySCj2yCSRww+i61YWe26ET47sydG9ONpSu2rShp4wEI/XkdBqSCpVdyByeyo3gxKytiT2zBIy9TLiZRNaAdDG4BV6QuwV3VtRu0exR6YbcdqIoLi/Pi8Jrx7yN76755QSA7Q4gLgNXbegx/yjk01PFE9vRvXkvt2WtRrgRvXgRvaA4ABdx6CcLZ1PQYrQDNjmj+oR0gGOXglFRD8WmR58XblyZC+o0Ny3lLgFXTeQ0GpIqI6PbQpV1wmj0vS7lcVfBiVkrSKStTpiRRuanLYdB3AYwWIvKLLH6ws52gSNqRTnh0jZytpHjvyJJcoQ5Fo+y3ZN0eRbS3LW9frltgReF/zirptGEz97Hiw2Gl1WLwFICJkMFHtkivHNX/DPCJOjem1iz4Yec8ReElXX7bPqijtv7Z6wQsvYJWStEuI1E3rRgipH8g4DkCN8BbQEn6jdY4MGmU5koSdm2gmhJ2Y/yt9Jf/WEqOuTunNhAooFKLbtEXZq874s9FqCsG+Cyuk6NdF6Int26w0232wvoSdvI0SC3bihoNgjU4ivclsUdfsHw7aldS0gbkGNG9BVZ4p/HHUkXHv31nLSus3UrlVFslFvCj0lqE5PRPgkwWcXALPsrjpgGc6ybWcRMmiafxIpgKo51tFa3FmeEjn/udOQVtcgmN01WudIL0W6bV+wVcBSAUsBTKWl0/zL6LDdBGD5o3n+B3QSdv79hJAjgWKPDJYgJ7W2tK4FJW5CjTUQ09pFnhB+SbRSuynXBi1RayBesqAU7WChJwu+cTSje0YZqJeBugHUTcCw3csKxR45CjRPBQWIKYCuAgkDUMUUoyD3Qbkr1y/2fBE2xW5F9toFnzeuJ2J9bYldS4VtqoChBgfiOum35n0bsC3fjl6Czy/+aG1DOsBu3FBQ7JHB4o9OyE5qruBTYjbUuAEt1oCuNlo+nfBH+KpIooqkXXG6dM2qY4FWsp3UrKjH8y9X7FkTgFUAzAJQrQHlKlAzgartZHI5n58MGtGLIf+tk1CAuOIEwjXFGU+puKNUEEO7LvKLPMCrh5q9FN5InbcNozlspdmyYUGDabv3LQ2WqQCGWOhvCREq3qNtwjsFuoHOSjGoaYNCj5CpgGKPDJZOjmqSobqiW9B0EzG9gZgr9mQLp7gr9pKouaNYqkjaVcTNOrSa2aq/E6uD8LMKQLUI1MpA0QBKFlCxnbI9jtkjRwMxoEjuU7LcqLJqAroBKAagGoDWFEzo/cX0NXrYqmOVZkliTvTkelfM46NhQIdh6TBMHZapedPG/a5mvaEtiT1hQu0Xfv7uElkE8mwkXWDNXigo9sjg6GSQ7ovwqboNLWZC1w3EVFnk1d3Z/3VPZC9l15wUrlGHVrVaA5LFOJUSWo0YRadGD25Er1YGipVWv0YFTkNu46j+jyHHIvKpkIAU/ILzn5gFGCagGYASFL2Tn8j/s1tiJ84vS2354spCz2wTe60xyw3oMGxJ7DVU2A0FqFlAxQJKNlBVgYYKNJQuYs9u1g22bBVlYSeLO79K9E+OJoRMBRR75OjiF4A6oGg2VM1yltpKLIm4Q8vZ03UCMBuI103oZRtaBV6xJws+d4lmjHLNiegJoTeBlmMaxR4ZJOIrLxpu5Z6Lnm0I/ui4NFvPs3TAjgG2Dpiaioaiu38uxZqx8rrnzyjZW8M5u+pWDPVGDI1aDGbDjeztPwz8chTYWwUyJwPZYaCmOCeOf9UBNGy3MUMIvKq0UxZ+8hJnoUcCE9IZztkLRUQ/FpnRiAtX00LXmfqvqa4vpxuH0JpCz/BelgwD8aqFWMWCIvvcyreS6DNLTkSvUnVSt0LoMbJHjgayRhOThoTQk6vVArO1QUJP1PIFDCu3Y4AVc8SeoehNEdfoJO782+w4GkYcjXocVl1zInv7DwM/exX47RhwFoCFJwF1tV3kifuGDVgifVuDV+x1Enryn1wUeqQP2KARCoo9cvTxmW1YoyUYrx2AmR4D/k8F2umWKBdvRvisio09r2fx6h4dC3PA2SdUkLGr7RecqrTc7WYDqDeAqunU6FXgXRR7ZBCIr7ioWgDa41YeK1wF0FRAlZszghxn2jylnWUnASOuwohrqOoJVJWkKHxo3jpDi+Ttzs8V21k1MwGjHodV1mGXVaCsAFYcSM8BchqgpZ2onji/KhWgdBgolIDqEGDMcbtw5ZOxBudPq3E4Z5xI1cr1eqzRI2SQUOyRacf8zQTs519DPf5n2FcloZ0el/sDocNArajgVzuOx6/+M4NLFys45ZIxnHRioVXu4w8WSMtsAHULqNneOIN8KaLYI1ONEHI2goMFPmtbxFQgpgNaDFCDxJ1f4KUApN2VAuyUAiOpoxKPo6InUVaFz4w8kjwlOU1L48rtNCpWGtVGCkY1DrscA0qKEx1PHAecsRDINAB1yBF7zb+WikDxd8DEHwBjIWCm3eLBGlp/SlXhCL0D7n2RyJYnmlPokZCwQSMUFHtk2rELddhvFGHpBWBcgQq9bUiE3bAx/qc43vhtBn+em0StpAHHoX3sg7/u23AySg27vR/QvwiZSsQ1Q7az9Y+ZjMMZvZJUgFgM0BOAFiTqZGGX8d63M4CVUWCkNVSTcZRjKZTUtOsgnZHcpH3LFvfTKJtpVBpp1KtJGOU4UNRaHe1KBjg+43wAMa+yWSNrAOUCUDkE5wF199PKdRUVd98EHIGXgNfnjVE9QgYNxR6ZdvSzhqDn34xkMgftbaLWp+XlCQCZXAPnvnM/Tk3tw+J5o5g7r+y9PvitmwRWyxmj0+WElxkylYgqheYcPbSCcm4QDhkAWQAZFcjEgGwMiGUcXYWstPIAcu7tEIA57q103xhSUc/pqGbiKMYyKCgZFJBrrgnkMeHejmPIXXlMiFs7j1Iti1oxjcZEEtaY5lgNHkZL3PmXGFZeygKNM903cxxaAk6eg1SCI+hiaLWiiHh6v7NlCPHBBo1QRPRjkdmEdlYeiQtPRyKXh64dBPAn94LZGgmbzjZw7vKDOP78P2FOsYQ5ExXH99bvzekXgO42jmglRws5PevPvnqEHoCM4gi9TApQ0vAKvRx6Cj3MAcycilo2hlI6iQLSrrDzr7wr8JzbcQw5ItDOo2DnUa1lUS1kYIzFW57S3QSf6HAqZ4HGGQBOR2vApV/oCbEXh3MGivStqGAkhAwaij0y/WgKlLgKJDS0DxFzUBRAj1lIxEzEahZUlbKNzCw8zRZwxF0CrSysLPSyKpDVgUwciGccoafIIk9E9YbcW5/Qs4cAY0iDkVNRTidRjKVRUNJuFE+O4Ik1p3l/wnaEXsHMo1zNoVrNon444Ub0FEfoiSULP09ED06WtmY5BtOowzv3qIhWGtc/eoUNGWQKYM1eKCj2SLSYxPUjWF4S0j/yNCF/6lakb0XgLgMgqwG5BJBKAVpGEnpy+lasDhG9Rk5DJRtHMZ5CQcu4EbwhKXLniLvDmNMUe4eliF7RzKFazqI6noE1psEaU7sLPU9ED46Oa5iAVUG7hY1Qg3IblOh7Z0SPkKMNxR6JFpNQbowvkCPBP+tYCL0U2oVeVnXq9NIJIJlyonoIiuh1SN1aQ4A1pMLMqaikEygmkijoGV+q1h/Rm4PDtiP+JpBHwcij1MihUsmiPp6CcSjeO6InhJ7QcTXTEXqmaMQQO0REzx/VEzV67LwlUwTn7IWCYo9EC15HyFGkWzOG3ECbhVufpwO5OJDIALq/GaOPGj1rSEUtp6OWjTmpW7W/iJ5co1eq51ApZlGbSMMc01vCrh+hJzScWXcjen6hV4J35IqI6MmpW56khBxtKPZItGBOlhxFJNe/wGYMT+pWBbJxIJtyZhPDn7rtIvTsIWeZeRW1bBzFjNyM0YroTbipWjmiJ4RewcqjaOVRqWZQm8igcSjZv9ArwPXGdYdVog6vwJObMXxTzVFDD0M4QsLDmr1QUOyRaMGaPXIU8A9E7jZeRTRjZN2InuJP3YrVoRlDjFdp5HQndaunXaGX90T05NTt4WY0z43oWXmUKzlUKhnUx5MwxXiVfpoxRESvZgGGsDcL6rr1R/T841UIIdMFxR6JFqzZIwNGtqoNSt36x6tkVSd1m0kBalAzhojodWnGMNxmjFIyiYLiTd1OSKlbR+jJYk8ar1LJoDqWhXlYhx2mGUME62oWYFfhHbHSKaJXhRP9Y0SPDAjO2QtFRD8WOWahciMDxD9eJagZwxPR05xmjEQG0OXUbR/NGLYb0TNzGsqppK/rNt8xondYHq9i5FGu51CpZtEYT7hCTwtXo1czgIYB2DW0p27lRgx5xArHqxAyk6DYI9GCOVkyIHqNVxEuZs2Inu6MV0m641W6NmPIyxV99pDipG6zcRRjKRS0bPhmDKPVjGGNaeEies3xKgZgyTP0eo1XEZ6FbMYgA4TduKGg2CPRgjV7ZEAEjVfp2HWrOgOT0ynHBq0totejGcMaUmDkNVTScZRSKanrts/xKmYeRSOHStkRekaYZowinGaMmuWmbmvof7yKcMag0CMDhg0aoaDYI9GCNXtkignldas49Xm5mDNDT+0U0fNH8yTRZ7rjVWSv204RPX8zhvC6LVdzqJayqE+kYHUbr+K3QBM6rgq3GaOK4K7bXs0YPKsImUlQ7JFowWsMmWLCeN1mVUfo5eW5K0GpW3/6tqfXbfB4lU7NGJVaFrWJjDswGeGaMSoAasLDtt9mDI5XIUcZRvZCQbFHogVzsmSKOBKvW6QBRU7ZdhqvIqVwG0OaY4GWTvTlddsUeQFet41+vW7bmjFMyes2SOgFRfQMcLwKITMbij0SLVizR6aAyXrdplMBqdtO41V8tXpNr9tEqm9nDK/XbR7VcubIvG7rZpfxKvS6JTMIjl4JRUQ/FjlmYc0eOUL6Ga8ie91mJa/bwGYMkbINEHrtXrcpn9etHM3L9/C6zaA+kYJxKAYc7lPoebxuDcAKqtGj1y0hsx2KPRIteN0hR4DcjBF3V8euW3i9bjuOV+nQiIE5gJ1XunjdytZn/jRuN6/bEKlbkZnt6nXbabwKu27J9GGrgD2A+jqbo1cImQUwJ0uOANnrVhZ6bV63SsvrNpPyDUzu1ozRnKHnzNEzchpqmTiK6aRbo+eP6M1x07buz3a+KfQKlpO6PXKvWxvhvG7FeBWmbgmZLVDskWjBmj0yCcKMV/F73ar9et3K9md5yes25ve69XbbymncCSmiV67kUCk7qdtJed1WbcAQAq7fZgxRo8doHpleTA0wB6BgzIh240Y0YEmOWVizRyZBt2YM/3iVjOt1m0s5Yq8vr1vfrZF3mzHSSRRiGRSaQi8fHNFrbnfHq1h5lCpZVA9n0TiUdGbp9YrojcOJ6BXhjlexAKPqbhRzV+SInkjfyuNV2HlLSL1ex4YNGzB//nyk02ksX74cTz75ZM/H/fjHP8batWvxlre8BZlMBmeeeSbWrVuH/fv3Bx7/7LPP4qKLLkImk8Hw8DBuuukmlEqlSb1nRvZItOA1iITgSLxuNb/XbY+Inj3kRPQMj9dtti11ezgooiePV6n5vW5DNGNMyutWiDxG9MjMwRpQZM/qI7J3ww034IknnsD69euxYMECjIyMYNWqVdi2bRsuuOCCjo/bsGEDxsbGsGbNGpx11ln43e9+h3vvvRdbt27Frl27MG/evOaxu3btwmWXXYa3ve1t2LRpE/7whz/grrvuwm9/+1ts3bo19Oei2CPRgjlZ0idH4nWrdxuv0qkZQ/a6jXu9boXYCxJ6nmaMhmjGyNDrlpBpYOfOnXj00Udxzz33YP369QCA66+/HosWLcItt9yC7du3d3zspk2bcNFFF3m2XXHFFVixYgXuu+8+3H777c3tn/vc53D88cfjqaeeQiaTAQCcfvrpuPHGG/Hkk0/isssuC/W+mcYl0YI1e6RPRERPFnpiWLKn69aN6GXiQMpN3WohXDGsIcDMK62ByckUCrEsJtScZ3Zem9Cz3Yie7UT0SvUcyuUsahMpGIcSbupW6bPr1gaqJlBvAGYNsEXXrVyrJ9RgBd4aPdkGjZCZgakpMDR1ypepdb8iPP7449B1HevWrWtuSyQSWLt2LZ577jns3bu342P9Qg8ALr74Yhx//PF46aWXmtsKhQKefPJJXH/99U2hBwAf/vCHkclk8Nhjj4X5XwWAYo9EDdbskR7Iqds4Ws4Ycn2ex+vWtT9LCpEn5uf5o3lzEFijZw2pqOZiKGaSKMbTKCi5tmaMwNStVKNXquVQGc+ifigFcyw2Oa/bRh2wi2gV7oVpxuBZQgjgpFcXLlyIbDbr2b5s2bLm/jCUSiUUi0XMnTu3ue1Xv/oVDMPA0qVLPcfGYjEsWbIEv/jFL0K/b6ZxCSHHFL28bj2pW8nrVgk7XsW9NfOy123QwGR/RE/qvLVzjtdtVXjdJrziLiiidxjenouKPF7FH82j1y2ZnZiaBlOf+niVqVlwyhaCGR0dxfDwcNv24eFh2LaNffv2hXq9TZs2odFo4AMf+IDnNRRF6fg63VLFnaDYI4QcE8giLyh128vrFmG8buc4HbeNnIZKKoGinnaFXl6qz+s0XiXvNmPkArxu0V3oyR23bV63ckRPLHrdktmJpWkwtakXe5amoJvYq1QqSCQSbduTyWRzf788/fTTuP3223HttddixYoVntcA0PF1wryGgGKPEBJ5jorXrW8Fe93mA8eryNsnkEfBzqNg5VAtZ1Edz8Ac02H76/PodUvIUSeVSqFWq7Vtr1arzf398PLLL+Oaa67B4sWL8cADD7S9BoCOr9Pva8hQ7JFowWAE8RHG6zbTr9etXKcnpW97e93mPU0ZXuuzltdtseGMV3EGJse8XbeyyKPXLTlGMaHCxJFNQH7iOw088R1vFG9ivPv3f3h4ODBVOzo6CgA45ZRTer7uG2+8gZUrV+K4447D1q1bPU0Y4jVs224+p/91+nkNPxR7JFqwtZZIdHPGCPK6zQmv27TUcRtUoxcg9OBaoNVyOmq5mJO69XjdisjdkBTN8zpjFGxX6BVaXrdtET1xP2hYsgjWGQ13vAq9bgnpxDUfjOGaD8Y82174uYnLlnZOky5ZsgTbtm1DsVj0NGns2LEDiqJgyZIlXV/z0KFDWLlyJQzDwLZt23DSSSe1HbNo0SLouo7nn38eq1evbm5vNBrYtWsXrr322n4/YhN245JowdErREKFk7rt2Yzhet1mXK/bZAbQQwg9ewiwhhQYeQ3VbBzFVAqFeNodr+KP6A3BP15l3M5jwsqjYORQrmZRnUijcSgFayzWGq/ir9Vrs0GzgbLlRPWMGmCX0FKBcq2ef7xKHbRBI7MNExqMAaxe0cLVq1fDMAzcf//9zW31eh0jIyNYvnw55s+fDwDYv38/XnnlFZim2TyuXC7jyiuvxOjoKL7//e/jjDPOCHyNfD6Pyy67DA8//LDHMWPLli0olUp4//vfH/r/FyN7JFpw9ApBOK/bjOpE9Dxet3L6No/gZgzZ/mxIeN3GUYxlfF633s5bWfDJEb1SxanRc1K3eu9GDNnrtoKW160tvG79A5PpdUvIkbJs2TKsWbMGGzduxIEDB5oOGnv27MHmzZubx916663YsmULXn/9dZx22mkAgOuuuw4//elPsXbtWuzevRu7d+9uHp/NZnH11Vc3f77jjjtw4YUX4pJLLsGNN96IN954A1//+tdxxRVX4PLLLw/9vin2SLTgNYsgnNdt1vW6zbrjVZRu41UChB7mAKbbjFFKpVBQ2mv0JqRonjNTr5XKLdjOLL1qNYPqeNat0euQug1qxhCTU2qWK/SC5ugFjVepg+NVyGzFggZzABKmn7PhoYcewm233YaHH34YY2NjWLx4MbZu3YoLL7yweYyiKFBVb/L0hRdegKIoePDBB/Hggw969p1++ukesfeOd7wDTz75JDZs2ICbb74ZuVwO69atw1e+8pVJfS6KPRItmJM9ppms1208yOu2RzOG7Ub0ml63MX/X7ZBvzXEHKA95vG5LtRwqNWe8ilOj16MZQ0T0RPldzQAMQ4ro9fK6rYNet4RMnng8jjvvvBN33nlnx2M2b97sifQBwGuvvRbqdS644AI888wzk3qPfij2SLRgzd4xy5R73XYRek4zBjp63fqbMZzOW3/3bR4lw+9120dETwxMbnrdmpNoxqDXLZndTEU3bvDzRjPSTbFHogVr9o5ZRDNGz65bFcgqLa/buF/o+V0xOo1XyauopOMoJVMoaMHjVVp1enNw2JZq9Mw8Sg3H67ZeSMEYi3f2uQ0SeiXbSdvWLXe8ShmtnG6/41X4zSfkWIFij0QLXr+OOUI1Y7het7mYm7r1Cz3RjNElomcNqc54lYzbjKF0Hq/SKaJXrjkRPTFHL9DnVh6vIurz5IieWXcHJvvn6AU1Y3C8CokWVh+ds5N7Xkb2CJn5MCd7zNHL69bTjKEA2RiQSwGq3+tWjup1qNFzmjEcr9tiJul23bY3YxyWI3qu360QegXL53Ub1HXbsRnDbk1N8Xjd9mrGoNctIccyFHskWrBm75hhsl63iYzTddvWcdtzvIrrdZtOoBhLdxmv0mrGGJebMaw8ypUcqpUMGuOS120/41WazRhWB6/bbuNV6HVLooc1oJo9C2bvg2YhFHskWrBm75hApG5V9DleZbJet9J4FSOnBnjdtqduZc/bw1LqtmjmUSlnUD2chXm4izNGV69bq4fXrVyjR69bQogDxR6JFlRukSfseJWMO14lIbxuxc4+x6uYrtdtOZ10um47eN12HK9iOM0YlUrGqdE7PBmvW8PtuqXXLSEAYECFMYDInhFRYzGKPRItmJONNHIzRtxd3bxuxXiVRMpnf9bneBVLeN1mY07qVgvfjFFqet1mYPrHq8hiz+91W0BLvxkGvW4JkbCgD2ioMtO4hMx8WLMXaVS0onqy0Av0ulVaXrdtc/R6jFexhwBbeN1m4iilkygo7c0YvcarFM08KpUMahNpNA4lvenabnV6RTjNGFVbcsbwj1cJSt2KockUeoSQFhR7JFqwZi+ShPa6jTldt4kMoHTyuu0S0TOGVNRzOqo9vW59zRju/gk7j3I1h0o5E+x1GxTRk71uywjwug1qxghK3YpmDH6zSXQZXIMG07iEzHx4fYskob1u3fEq6OZ120HoCa/bal9et3PavG4nbMfvtlrNojqegTEW8w5M7hbR83jd2pP0uo1mGooQMnko9ki0YE42UnRrxug0XiUddwYmq528bjuMV5ms161nvIrpRvRcr1trTAPkZowgodc2XsUEjEYIr1sxXoVet+TYYXB2aYzsETLzYc1eZOjlddupGaOn122X8SqNbHivW3m8SsnIoVJqed1avYReoNetaMYoSatbM4YYr8LULSEkGIo9Ei1YsxcZgrxu5fEqzVgPcvgAACAASURBVMCd6jRkBHrd9jFepel1m1NRSSdQSiZ9XrdyNC8f3IxhOEKvXHYs0IxDceBwn6nb5ngVE6ibktdt2PEq/CaTYwcT2kBGrwwiWjgToNgj0YLXu1lPt2aMtvEqwv6sk9dtHzV6Ta/bbBxFPe3zupWtz/xpXMnrtu543dYm0rC6NWMECT163RJCBgzFHokWzMnOekTqtpPXrSd1q06V12080OvW74jRjO4Jr1vLacbo6HXbsxnDbuk3et0S0jcWtAHN2WNkj5CZD2v2Zi1H4nWLNFop2ynxum1F9PxCr2+v234iejULMOh1SwgZLBR7JFqwZm9W0svrtq0Zw+d1q3SK6PXldZsM9LoNjOjJXreW3+u2j2aMUF63QRE9et0SAji1dYPpxmVkj5CZD5XbrKMfr1u5GSPr97oN0YwhxquYOdUZr9Lmddvqtu0Y0TPyKNVzqFQzqE8k271u+2rG6MfrVhZ6dTgRPXrdEkLCQ7FHogVzsrOKKfW67asZQ3GcMbLxQK/bTqnbNq/bYh9et50iepP2umXXLSECOmiEg2KPRAvW7M0qxHiVvrxuVadGL90potejGcOSvW5TSRTUPpox7DwG53XbKaJHr1tCyNRCsUeiBWv2ZgVH4nWrdhqvMhSw/F63GdfrVgmq0Wu/7/e6rZYzqI/36XXrT902vW77Ha9Cr1tCOkEHjXBQ7JFowevhrKBbM0ab163S8rpV/ONV5NStiOoFet2qrtdtsqvXrYji+VO3Hq/bQ/GQ41Xg6LiaDUfA9duMISJ6bMYgxA+HKoeDYo9EC+ZkZzST8brNyF63Icar2O54FUM0YxyB1221lkVdeN2GccYQXremGK9Cr1tCyNGHYo9EC9bszVgm63WbSgU4Y/hr9PzRPOF1m9OOyOu2aORQLTsRvb69bv3jVRoGYFUQnLrt5nXLiB4hneBQ5XBQ7JFowZq9GUsor1sVSMeBZD9et20dt16v22Ii6Ruv0hqcLCJ6nbxuK+Us6uPT5XVLCCFTA8UeiRa8Rs44jsTrVu+3GUOqz7PyktdtLB04MDkookevW0JmD2zQCAfFHokWzMnOOCbjdZtNAVq/zRiu0LPdn7t53YpmjMNyRA9DTa/bguWsajWD6kQGRpjxKh6vWxv9e92KRa9bQshgoNgj0YI1ezOGI/G6Vfyp207NGLL9WV7yutX9Xrf+hoyAZgwrj5LrdVsfT/bndTuBDl63nbpue41XIYT0gzUguzTW7BEyG2DN3owgtNet6jRjZMJ43fpum163yWTbHD15vIo8OFluxihYeVQrGVTH6HVLCIkWFHskWlC5TTuT9rpNA7qcupWbMTpE9OwhwMirMHMaymnX61Zrb8Y47IvoHZYjeq7XbbXqRPQcoadNwuvWAKxuzhhB41XYjEHIZDChDmjOHmv2CJn5MCc7rXTzuk3CieIFjVdJBo1X8XvdBjZjuF63ubiTutXCN2OUDKcZozrhjFcJHdGrwhV69LolhMxMKPZItGDN3rTSzeu2res2jNdtQDOGlZe8bpP+gcl9jFcxnTl6lbLTdRuqGUN43dYspxkj0OtWCD163RIy1ZgDmrNHBw1CZgOs2ZsWJuN1K8ar9O1162nGCPK6zfuszgKaMdq8bp05el3HqwR53Zbhet3K41XEkpsxgubo0euWEHJ0odgj0YLXz2lBrtMbiNetr1bPzKmo5lyvW3gt0CY8dXpOM8ZhX+q2YOdRrR2p1y1Ar1tCpgd244aDYo9EC+Zkjyr+ZowkBuB162nGcL1u00mnRq8topdHX1631UF53XYar0KvW0LI9EGxR6IFa/aOGqG9brUj8Lp1b71et96Iniz0OnrdmvS6JSQK0EEjHBR75JiHsZbJ0a/XbUZ1onrpREivW3dZQ44FWnevW//QZMcVQ6RuC0YepUYOlWoW9YkUjLE4MKZ4RZ686HVLyIzGhDag0StM4xIy82GYbuCE9brN6UAu3sPrVgg9EcmT7ltDwus21qfX7VAzqjcBpz6v5Pe6DYrmyc0YIppHr1tCSASg2CPRgmncgaOiezOGJ3WrANl4H163AULPHgLsOeG9bkVEb1zyuq3UMqhNZNDwj1fx3w/0urXdZgx63RIyU7AGNHqFDRqEzAY4emVgdPK6DRqv0uZ16xd6eQQ3Y8j2Z0MqGlk9tNdtM6Lnet1WXK9b0+9126k+r6vXbdB4FXrdEkJmNhR7JFrw+joQennd+serZFQndRvodesfrxIg9DAHMNxmjFKfXrfyeJWC7TRkVCsZVA9nYY75vG57RfRE6rYW1uu2Dkb0CBk8bNAIB8UeiRbMyU45/XjdeiJ6stdtRtoZ1IwRkLo1hlyv21SQ1+1Q22rzujXzKNVyqNZERM/ndStH9ITQ80f0qq7XbaAzhtyIIYReHRyvQgiZqVDskWjBmr0ppdd4la5et/7RKj2EnjNLz/W6zcZRjKVQ0LKeiF6rNq9d9DW9bhtdvG6DmjHG4UTz2sar0OuWkJkKhyqHg2KPRAvW7E0pQeNVQnvdyjP0ukT0rCHX6zYdRynl97r1pmxb41VESlfyuq1kUSsEeN12S93KXrd1C7BqcESdsMzoNl6FXreEkJkNxR6JFrzWTgmhvG4Vpz6vo9etaMboEtEz3fEqXq/bHCakVG2b0JOaNZpet6Us6oUUrEN6Z6HXqRmjCp/XbT/OGPS6JWQ6sAZUs2exZo+QWQBzslNCKK9btU+v2w5CD83xKjGU0u3jVQIjerIPrut1W6llUZvIuAOTEWK8CtxmDBv0uiWERJFoSlhy7MKavSPCX58ne92KRgy5GSPnztCLux23YrWlbTukbhtDGsq5OIrpFAqxLCaUHCaUPCaUfHNu3uGgiJ495Ig8M49yOYfqWBaNMdfrdkxpuWN0G68iyu9qJmBW0Jqi3KsZg163hEw3huugMYjVi3q9jg0bNmD+/PlIp9NYvnw5nnzyyZ6P279/P2699Va8613vQj6fh6qqePrppzse32g08JWvfAVvfetbkUqlcPLJJ+Oqq67Cvn37Qv2/AhjZI1GDNXuTZmBet13GqzS9bhOpQGeMIKEnR/SKZh7VcibY67bf8Sp1s8t4FXrdEkK83HDDDXjiiSewfv16LFiwACMjI1i1ahW2bduGCy64oOPjXnnlFdx1110466yzsHjxYjz33HMdjzUMA6tWrcKOHTuwbt06LF68GGNjY/if//kfjI+P45RTTgn1nin2SLSgcpsUIm07EK9bn9CzXAu0ltdtyud1OxQs9ORmDOF1W8mgPp6CcSgGHA4Qer28bhuG63UbNF6FXreEzFSmy0Fj586dePTRR3HPPfdg/fr1AIDrr78eixYtwi233ILt27d3fOx5552HP//5z5gzZw7+/d//vavY+/rXv45nnnkGP/nJT7B06dLJfRgJpnFJtGBONjRyM4acupXr8+QavZwO5JNAOu2OVxECr4+uWzFepZbTUcomUEykUVCzkjPGHHduXnBETzRjFOs5VApZ1A6lYYzFgMNK74ie3+u2UQeskrSj34gemzEIOVZ5/PHHoes61q1b19yWSCSwdu1aPPfcc9i7d2/Hx2YyGcyZM6fna9i2jW9+85u45pprsHTpUpimiUqlckTvm2KPRAvW7IVGRPR0AHH08LpVHa/bTApIZgA9pNetNQQYec31uk2hEE9jQu09XmXcblmgFY08ylVnjl7jUArWWKxVo+dvzGgbmmwDZcuJ6hl1wJbFntyM0bTQQGu8iiz2CCHTiXDQmPrVXRbt2rULCxcuRDab9WxftmxZc/+R8utf/xr79u3D2WefjRtvvBGZTAaZTAbnnHMOtm3bNqnnDC32SqUSvvjFL+LKK6/ECSecAFVVsWXLlsBjX375ZbznPe9BLpfDCSecgA9/+MM4ePBg4LH//M//jLe97W1IpVJYuHAh7rvvvrBvjRDW7IVA1OcFRfT8zRjC/iyXdrxuVTma50/bdojoGUMqqrk4iukkCrFW2rZXRE+kbids1+v2cAb1sQCv205pW894FRswamhNUe7X65Yij5CZhDUQoaf1TOOOjo5ieHi4bfvw8DBs255U84Sf3/zmNwCcVO7TTz+NBx54ACMjI6jVarjyyivx4osvhn7O0AnvgwcP4ktf+hJOP/10LFmypKPK3Lt3Ly6++GIcd9xx+NrXvoZCoYC77roLL774Inbu3Aldb730P/3TP+FjH/sY1qxZg09/+tN45pln8MlPfhKVSgWf/exnQ38ocgzD63HfhPG6lTtvlfSRet2m2rxu+xqvYklet4djsP3RvE6dt0LTVUGvW0LIEVGpVJBIJNq2J5PJ5v4jpVgsNm9feOGFZjPGpZdeigULFuAf/uEfOgbZOhFa7J1yyinYv38/5s2bh5/97Gc4//zzA4+74447UKlUsGvXLsyfPx8AcP755+Pyyy/HyMgIPvrRjwIAqtUq/u///b/4q7/6Kzz66KMAgLVr18I0TXzpS1/CjTfeiKGhobBvkxyrHOs52T44Eq9bLWQzRr9etx3Hq7jOGOVaDpVqFo3xBMzDutcCLUjo+SN6NQMw6HVLSFQwofY1JmUyz9uNVCqFWq3Wtr1arTb3HyniOS688EJP1+2pp56Kiy66CM8++2zo5wwt9mKxGObNm9fzuCeeeAJXXXVVU+gBwLvf/W4sXLgQjz32WFPs/fd//zcOHTqEj3/8457H//3f/z0eeeQRbN26Fdddd13Yt0mOVViz15Uj8brV/aNV+vS6beR0Z7xKvN3rtp/xKsLrthbkddspotfmdWvS65YQ4uGX33kJv/zOy55t1fF2ISczPDwcmKodHR0FgNAjUYIQz3HSSSe17Zs3b96k6gIHMnpl3759+OMf/4jzzjuvbd+yZcvwgx/8oPnzL37xCwBoay1eunQpVFXFL37xC4o90j+s2etKJ69b0YgR5HWb8o9X6aPr1nKFnpHXUEm7qVtPRM+pxes6XsXMo2TkUC5nUZ9IwTiUCK7R6yT0PF63VfTvdSuPVzmWvh2EzB7MKRi98vYPno23f/Bsz7Z9P9+P+5f+S8fHiPK1YrHoadLYsWMHFEXBkiVLjug9AcDZZ5+NWCwW2Nm7b98+nHjiiaGfcyDduELhdipiPHToEBqNRvNYTdMwd+5cz3GxWAwnnHDClBQ7kmMIXpsDkVO3cbQ3Y8ipW+F1m3c7btvGq+TRM6JnDamo5mIoZZIoxjOOM0azGWPIF8WTxqpI41VK1Rwq41nUD6VgjsU6C71xBA9MrsIZr2IX0d94FX8zBr9MhBAvq1evhmEYuP/++5vb6vU6RkZGsHz58mY2c//+/XjllVdgmmbo18hms1i1ahWeffZZvPrqq83tL730Ep599lmsXLky9HMOJLInChR7FTHGYjFUKhXE4/HA50kmk1NS7EiOIY6lnGwIJuN1m5fnrkzG6zYXQzHlet0q3oieE9Ubat7v6nV7KN5fRC/Q67aO/psxamAzBiGzA9GNO4jn7cayZcuwZs0abNy4EQcOHGg6aOzZswebN29uHnfrrbdiy5YteP3113Haaac1t3/5y1+GoijYvXs3bNvGli1b8MwzzwAAPv/5zzeP+8pXvoIf/ehHuPTSS/HJT34Stm3j3nvvxdy5c7Fx48bQn2sgYk8UF/ZTxJhKpVCv1wOfp1qt9lHs+F9wLlkyi4D/v717D5KiPPcH/u2ey851F9mVw0KEHIObeMMVIiGSCwdNQDSkjLAeKwG0CP5KJRhjBE0kprwFA554FEVNLFdEoxE1MRJNHT0SwgFERJJ4JRpFAwsGYe9z6+7390dP93T39MxOLzvszvD9VE0t9nT39jJZffK87/M8ONXlXKp63LNnYw3ynJMxjJetGMOvL90Gs5v3JCOjF7e86lCw6jZT59NHoEVq0O2PoEuKWDJ6tY6snt5upQN19mKMZBzJZAyZ9uys2/YSq27NYgwVUNPQAz1nRq/YrFsNzOYROf0NgLPVR3IwHmTIePjhh7Fs2TKsXbsWhw4dwvjx47F+/XpMmTLFPEeSJMhy/uLpT37yE0iSZJ5jBIiSJNmCvRNPPBEbN27E0qVLccstt0CWZZx11ln4+c9/7rpq2peyBHvGgxjLuVZtbW0YPnw4AoGAea6qqjhw4IBtKTeTyeCTTz4pYbPjDADef3CqUtyzZ+rvrNtIGJDLNOvWWMbV9+kNy/7ZMutWq0UiO+tWPeT31l6Fs26JyuBU5CdP2gDc73LukaNlmyqX4759CQaDuO2223DbbbcVPOfBBx+0ZfrM+2ul/3umubkZf/zjH0s+v5iy7NkbNWoUjj32WGzfvj3vvW3bttk2MDY3N0MIkXfuK6+8Ak3TBmSzIx1FqjVy88jaQ88to2cN9GKyPgLNmHUbyO7Tk9x66BUoxlDqZKTi/uys24g+61a2L90agV1HNtBrF3X6ZAyRnXWbjCPRlZ11eygAccgHHCow79Zt1m1SBTIpQEtkJ2MUmnVrnYxhFGMw0COi6lW2cWkXXHABnn32WVs1yYsvvohdu3ahpaXFPDZt2jQMHz4cq1evtl2/evVqRKNRnHvuueV6RKpG1bwmWyLrrFtnMYbbrNuYH4iHgYhbMYaXWbdxt1m39pd7P71adGccs25LzehZB18UnXXrtkePxRhElWqwxqVVqn4t4959991ob283A7lnnnkGH330EQBg8eLFiMfj+NGPfoR169Zh6tSpuPLKK9HV1YWVK1fitNNOw8UXX2zeKxQK4aabbsKiRYvQ0tKC6dOnY+PGjXj00Udx6623ljQ0mMjEPXuQYa+8LTjr1tJeJRoG/B6LMYSlvUoyGkRPOJTdo5ef0etwZPSs7VW6tVokElGkOqPIHAx5LMYQ+gg0sxjDmdFzC/SMpslcuiWio0O/gr2VK1fiww8/BKBvKnz66afx9NNPAwDmzp2LeDyOT33qU/jTn/6EH/zgB7juuusQDAZx3nnnYeXKleZ+PcNll12GYDCI22+/Hb///e9x3HHH4Y477sD3vve9w/zx6KhzFO/Zs2b0nHv0InApxgjor5qoY49eicUYSp2sN0yOBNEdiOpVt67FGMMsS7h1tvYqvck4Er1RpDvDhWfdurVXyZt1m4L7Hj3OuiWqRip8ZZqgMfD3HAr6Fey9//77JZ134okn2hooF7NgwQIsWLCgP49DlHMU//e7WDGGs71KNNteJRYGpKhj1m2JxRiqMes2bJ11a2+vYl++NYK9XHuVZDKGZHuUs26JiMqoLNW4RIOm2tZkS1Bs1q1rexVj1m0U8FmXbkvI6BmzbhVj1m3AWnVba2bx7Pv0shW3brNu20ucdZvXXoWzbomOZtoATNAodN9qxGCPqstRtmfPc3sV66zbYu1V6hwvsxgDJc+67XT01TNn3SrOWbf9aK9izrrtsbw465boaKGWqfUKCzSIKsFRtmfPbdZtGC6BngzEpAKzbo3qW2vVbd74M30EmlorF5h1W6Da1jnrNtOPWbdm1a3brNtC7VU465aIyMBgj6rLUfLf8mLFGM4+elFJ358XD+hBns+Z0bMu3boEesas21Tcj1Q0W4whWRsm5/rnOduqWDN6vSk9o6cXYxSZdVto6TYBfTKGSMK96ta5dGtk9NhehajaDNa4tErFYI+qSyWvyXrQ16xb29JtthgjHgbkYu1VrFk9y7KtOes2FkB3NDvr1lKM0VEo0BOWYgytFomkMevWQ0bPaK/CWbdERP3GYI+qS5Xv2TucWbdSBLkl22LFGJavSq1l1m0g4mivUjijZxZjaLXoTcSRTESR6TBm3cJjMYZWZNZtofYqnHVLVM1UyGVqvcI9e0RDXxXv2TOWbo0xaGWZdev42tes276WbrvV7Kzb9li26nagZ90Wm4zBjB4REcBgj6pNpURuHhVrr+IM9GztVbKzbl2LMQpk9LQ6QK2VocZlPaMXDOuzbl0mY7Sjzj2jp+jFGIlEtmFye8Bbe5Ue6O1VMmq2GKPQrFu3hslGMQYRVSsV/rK0XinHPYeC6vyp6OhVSWuyJXLOurWOQMsrxkCuvUpNGPC7ZfPc2qtYAj5RKyEd9yMVD6LbH7FU3Vozd3W2r86MXo8x67YzCrU/7VWS0PvoaUZQ58zoFWqvwmIMIiInBntUXapwz16/Z90WqrotEOiJOj3Qs8+6dc/o2V6O9irdqjHrNnIYs261Ag2TOeuWiACtTH32NO7ZI6oAVbRn73Bm3Uqlzrq1FWPISMf9SBrtVRDNFmNYq25z8277nnXr7/+sW+F11i2LMYiICmGwR9Wliv57359Zt/Fsqq+kWbeOr2rch2Q82zC56KzbYcVn3XZEoRwKAF5n3SYApEU20HOrui0261Yd8L9/Ihq6OEHDGwZ7VF2G+ppsCTzPuvUDkWx7FdljexVRp2f0lLgPvZGQvkdPLjGjZ511m4wjkdJn3WqHfIDnWbcqoGYKZPTcZt0a7VU465aIqC8M9qi6VPievf7Oug2HXSZjFJp1ay3GyJt16xyBZq28HYaObJuVdues257crFut1EDPNutWscy6LaUYw2ivwmIMoqMRJ2h4w2CPjnpDKVRwm3XrWnWbnXUb8Trr1hx/BmjW9iqhkCPQq0VnXlavDh0it3RrzrpNxJDuCkM5FLQv3bY7XkagZ67QGrNuVc66JSIqIwZ7VF2GUprOAy/FGCXPujUCPcuM21wvveys21i2vYpUuL1Kbul2WC7QE7XoMWbddoWhHgy4Z/OsxRhGNs+a0eOsWyLqB07Q8IbBHlWXCl3GNZZuC826tRZjxCS96rbPWbcugZ591m3QddatUYxhG4Umcv30ukQtOrVaJFMxJDsjUJztVZx/dm2vAj1+46xbIuoHFb4yNVXmMi7R0FdhrVcKzbot2F4lO+u2xph1W0p7FUvQp9QVm3VbuBjDzOhptejJzrpNu826LbQ/rxO57XcpDVA465aI6EhhsEfVpYJigWKzbp1Vt1EA0f7MunX20ovL2Vm3IddZt27FGNaMXpemN0xOtsegHnLMuu0ro2fOutU465aIDgsLNLxhsEfVZSisyZaglFm3toyec9at8aZbMYbL0q1Spxdj9IZDRWfdGhk9s+JWZAM9pRY96TgSqSjSHSGohxyzbq0ZvQGbdZsGZ90SER0+BntUXSpgz16xYowQ9Diu5Fm3fQR6ejGGPus2GQvqS7dFZt06M3pGMUZ3Ri/GcJ1161aMYVTddiG3KstZt0Q0QNhU2RsGe1RdKmDPntFexa0Yw9ZeRdKLMaJBfek24FZ120dGT6uzzLqNhNAl5Rdj2DN6w9Au6myzbrvUOBLJmPus22JLt0Z7Fdus217kRmZYAz1nexVj1i0DPSKiw8Vgj6rLEI4LPLVXyY4/M2bdytalWyOb10d7FaXObdZtvM9iDOes22RvzH3WbaH2Kq6zbou1V3Hro2cUYwzhD5SIBo0KX5lar3DPHtHQN4T37BUqxijUXsWYdSt5bK9ifFXjMpKxIHrCoYGfddtnexVkV2UF9ACu1GIMzrolIhpoDPaougzBPXuHPevWbTJGgfYqRjGGEvfpxRiBsKXq1mUiRpFZt8lUDGm3WbcltVdRASUD90CPs26J6PBoZeqzx2pcokowxPbsDfis21r02V4lE7POuo2VVIxhnXXbrcSR7Ikh6TbrtlBGzzkZwzbr1niVMuuW7VWIiAYagz2qLkMsKeQ269baXsU261bWM3qhUmbdOgI9LTsCrfis2z6KMZRa9ChxJLJ79JSDQaDdZem2aHsV9TBn3RIR9Y3VuN4w2KPqMkT27BUrxshrmCwBcT8Qzy7d+ss869atvUqnqEVv2mivEoFWajGGMQCDs26JiIYsBntUXYbInj3Ps26DQCwM+Mow69a5R8/cpydyI9C6tFokk1EkO6OHMetWoPRZt8aLs26JyDtO0PCGwR5Vl0Hes1do1m2xYgzXWbfW9ipFZ93K+h4911m3ddnii2HmMm5eMYZt1m0of9ZtoWVba3uVfs26tbZXISKicmKwR9VlEGOHYrNu3YoxorJejBENAVIUkApl9IrOuvXpxRihkOvSrdFWxViudRZjdGm1SCaiSB6KQW3320egFau85axbIhpEKuQy9dnjnj2ioW+Q9uyVMuvWTNw5Zt36B2LWrS+/vUq7JdDrMBslW4ox0nEkk9lZt+1+iEO+4oFeXsNkJVt1m4L7rFtnRi8NfY8eizGI6PCoZWq9wqbKRJVgEPbsWYsxgtnXkZt1m99epd2R0euAtUhDz+j1ZGfdGu1VSsroFWyv4mXWrRHoMdgjIjpSGOxRdRmEPXvWWbfOQC9v1q08gLNuw4UaJjsCPWFMzNAbJnercSQSetVtXjFGsaVbY9ZtStOLMVxn3RZqr8JZt0Q0cFig4Q2DPaouRzCO6M+s23ggOxljIGbd2vbo1RXM6LnOuu1wmXVrDfQ6kB/ombNuvbRX4axbIqLBxmCPqssR3LM3aLNuI4XaqxRfurXNuj0YLD2jZyTuEtBjN8+zblmMQUQDSytTU2WNBRpEFeAI7Nnrz6zbaKFZt8ar6KxbHxSjGCMQQZdktFdxm3XrUoxhzLpNWmbdWidj9LV0a8y6VY32Kpx1S0RUSRjsUXUp8549z7NufX3Mui2hvUrGaK9iq7p1TsaoLZjR61b1pdtkh8us21Lbq2RUQEuAs26JaChQy7Rnj9W4RJWgzEmkUmfdxmR9n16kZoBm3daE0OXva9atoxhDyVbdJvQ9enmzbvssxoBl1m0CnHVLRFSZGOxRdSnTnr2+Zt3mVd1mizFqDnvWbUBfupW9z7rtscy6LVqMUWjplrNuiWiIYlNlbxjsUXUp0549Gbl9em7FGLalW1kP9vo761b0Meu2Mxvo5UahDUO7yAV6XVoturVaJArNuu2zGIOzbomIqgmDPaouA7xnr9CsW7f2KkYxRsyYdetcuu2jGCNv1q3fWhEFvgAAIABJREFUOevWWZAxLDfr1sjoabXoTcSR6GvWbdFiDGPWrVF169yjx1m3RDS4OEHDm+rMV9LRawBjjWKzbp1Vt0ZGL16jZ/RqnLNuSyzGUGp9SMSz7VUCUVug12lZujXarLSbAV9u1m1PIobkoRgyB0PQDvm9t1dJaoCSRG5khjWzZy3ISFm+KmCwR0RHg3Q6jaVLl2L06NGIRCKYPHkyXnjhhT6v27dvH6699lpMmzYNtbW1kGUZGzduzDsvkUjg7rvvxvTp0zFq1CjU1tZiwoQJuPfee6Fp/Vs5YbBH1WWA9uxZgzxjKoa1vUoE9oxePABEQkBNRB+BJseywZ5z2bbA0m2mTkYyHsjOuo2g0xdDp2zvo9dua5ys/7lT1OlBnqJn9JKdzlm3UukZvaQCZJKA6IW3WbesuiWiI8uYoDHQr1ImaMyfPx933HEH5s6dizvvvBN+vx8zZ87E5s2bi173zjvvYMWKFdi7dy/Gjx8PSXL/D9Y//vEPLF68GABw9dVX4/bbb8fxxx+Pyy+/HAsWLPD+lwUu41K1GYA9e321Vyk06zbk1l7FLdhzBH2ij1m3JRVjKP2YdZvXXsWYdVtqexXOuiWio8u2bdvw+OOP4/bbb8dVV10FAJg7dy5OOeUULFmyBJs2bSp47ec//3l88sknGDZsGJ588kls2bLF9byRI0fi9ddfx4knnmgeW7hwIRYsWIDW1lYsW7YMxx9/vKfnZmaPqssA7NkrtEfPterWl5t1WxPVs3p5e/QKBHqiDlDrJGSss26DUVtGzzkKzSzGMDJ6ai2603H09saQzM661Q4FSs/o9QggqQKpDKCmslm9btiXbo1o0BihYWT1uE+PiAaHMUFj4DN7xcOidevWwe/3Y+HCheaxmpoaLFiwAFu2bMGePXsKXhuNRjFs2LA+f7b6+npboGc4//zzAQBvvfVWn/dwYmaPqsthxB2eZt1KQDzomHXrltErkM1DHaBm26u4z7p1y+jVwdpE2Zx12xNDujOs78/zPOsWnHVLRBVHgQxfGYoplD6CvZ07d6KpqQmxWMx2fNKkSeb7o0ePHvDnAoC2tjYAQENDg+drGexRdTmMPXvWyts+Z91m9+nVWvuuuC3dFtmrp7dXCbjOunVtr+JYuu0StUikYkh1Hs6sWwHOuiUiKk1bWxsaGxvzjjc2NkIIgb1795bl+2YyGdxxxx04/vjjccYZZ3i+nsEeVRePiSYjwPMhN+s2hNJn3SICSMZEjGLtVYz9edmK20zcp7dX6XPWrSXIc5l1mynLrNtCGT3OuiWioUGDvyytV7Q+7plIJFBTU5N3PBQKme+XwxVXXIG3334bf/jDHyDL3nfgMdij6tKPzF7ZZ906vpqzbmvCJU/GaLfNuq1Fsjd6eLNu0yogjFm3zkDPOQKNs26JiAAgHA4jlUrlHU8mk+b7A23FihX41a9+hVtuuQXTp0/v1z0Y7FF18ZB0MoI845dgQGbdFsjoaXWAVmuddRsuYdatXoyRP+s2mp11GwDaSwz0rLNuMwqgJcFZt0RUqYwCjcOR+PUzSPz697ZjoqOr6DWNjY2uS7XGfrpRo0Yd1jM5tba24tprr8Xll1+O6667rt/3YbBH1aXEzJ51+db4c8GqW2Tbq2QnYxRsr1KoGGOYHuil4n6k4gF9MsaAzLr1uHRrzLrVEsgP9Iq1V+GsWyKqPuGLZiF80SzbscyO13Fg4qwCVwDNzc3YsGEDuru7bUUaW7duhSRJaG5uHrDn+93vfoeFCxdi9uzZWLVq1WHdi61XqLqUEI9YAz23pVvbrNtse5VYUJ+MEYpY2quUUIwhslk9tVbW++hFwugKRkpqr9IuLLNulVr0JvU+epmD4dLbq3QB6BZAr6Zn9ZQ0IHrgPhmjUHsV7tMjoqFFLVPrFbWPsGj27NlQFAX333+/eSydTqO1tRWTJ082K3H37duHd955B6qq9uvn27hxIy666CJMnToVa9eu7dc9rJjZo+rSR2bPWXFr/Br6MICzbm3jz2Rk4v5cMYbXWbciO+u2N4p0Zwhqf2bdJp2zbkttr8Igj4jIatKkSZgzZw6uu+467N+/H+PGjUNrayt2796NBx980Dzv2muvxZo1a/DBBx9gzJgx5vGbb74ZkiThjTfegBACa9aswZ///GcAwI9//GMAwIcffohZs2ZBlmV861vfwm9+8xvbM4wfPx6nnnqqp+dmsEfVpUhsYgR6xhg0gVy5gQ/57VWisr50GwvrgV5Js24dXxWjGCMUyuuj5z7r1tFeRatFMhFFsj0GtT0AcajEpVujvUoSensVkYS39ioaWIxBREOVpvmgagPfZ08r4Z4PP/wwli1bhrVr1+LQoUMYP3481q9fjylTppjnSJLkWjX7k5/8xByTJkmSGSBKkmQGe++//z66uvS9g4sWLcq7xw033MBgj45yLpk9t2VbI+gTyDVStmX0fHoxRk0U8Fn76JXSXqVOz+gpcR96IyF0B40RaLXoRNylvUpu7q2tvUoqjkQyhkxHTXbWrYdijF4AKUUvxhAplDbrlu1ViIj6EgwGcdttt+G2224reM6DDz5oy/QZNK3v/xP91a9+td/Lv4Uw2KPq4hKjWFurBKGHM0YVrkBuWddt1q2/1Fm3loBP1ElIx/xIxp2zbmstwd2wbDsVe9BnFmNkjGIMzrolInJSVRlQBj6zp6rVWcrAYI+qizOzJwGyBPilXLCnIZe/MrJ+AQAxSV+6jQb1Pnqu7VWKBHpaHSBqJShxH5KRGvSEQuiSY5alW2ugl6u27RC15h69LrUWPYo+6zbdGYZysMbbHr0eAaQ0IK0BWgqFq26LtVdhoEdEVE0Y7FF1sqzdyjIQkICgBGhCD2UC0IM9s1hD0oO8WHbWbcH2KkagZx19ln1ptlm3EXQhagv0ctW21myeJdATtehJ6Rm9dFcY6qGAPdBrR/6s204M0KxbtlchosqhKj5AGfgQRi1DtnAoYLBH1cUaq2TXb2Uf4Jf1gE5kzzHKD/zQA7+QrAd6sb5m3brMuLXNuo0H0BPOzrqVcu1VrEUY1kAvb9ZtMjvr9lDQPZtX0qzbQiPQ3IoxUmAxBhFRdWOwR9VFgr2/ih/wBYCaGkCSAFkF/Fouh+WTgKAMBAJAwJh1W0oxhnX8WZ1l1q3fOuu2yNItatEp6tBlm3Ubzc26LaWHXidyq7K2WbdGD71SijE0MJtHRJVGU31l2bOnqczsEVUGR/mtLwjImh7Y+RUgqObCG58MBPyArwaQS51168jo9TXr1i3Q67As3XZpcSR7Y0h2xPJn3faV0bPNui3UXsWtGIOzbomocqmqDFGWYI8FGkRDkzWbZ5TZBmBWZMjZyE7yAZIC+JTcpXI28+ergfus2wJLt1p2j17hWbfFM3rGrNvujN5eJd0ZhnIoALgFesZevT5n3bq1V+GsWyKiox2DPaoO1qVbo+w2BD2myQaDkl9fxpUs7YskPyAb53rYoyfqpCKzbnPFGHkZPWPpVmQDvS7LrFu3pVtrMUYHciu0ZjFGJttehbNuiejooSo+aJmBz+yVI1s4FDDYo8pmtFpxZvVqkGsbl+2gLAUAn3Pl0pibZgR7fQR6IhvoKbU+JKNBvRhDiuRl9PICPWt7FU3fp5dIRpHqjCBzMJy/XFu0GEMASZGdjOHWMNmtGMOYdculWyKiow2DPapcxvKtc+nWCPSMuMaXfWVgr0eQYM8C1qLPYgylzph1G0R3IOqYdWttr+JYujWXcPVZt0lz1q3fPdBzK8boRnZFVgBKKhvoWYM7zroloqOD0HwQahlCmDKMYBsKGOxRZTMyekbAF0Qu2LMOvjVGaBirl0agaAR7YZRUjKFmizF6QmF91q1kL8bos72KVotkMopkRwzqIces21LaqyShN00WKeRX3XLWLRER5WOwR5XHrRjDCPKMlUprEss4z7pVTc7eJ5h9RVB06Vapy866DYeyI9CsxRjW5VtHw2S3WbftNVAPOWbduhVj5LVXUQBFsWT0+mqvkgZn3RJRVVLKMy4NCqtxiYYGZ7BnzegZq5Wq5VwjIFRd7mFcF0GRYgzoS7exILqD1lm31sxdnSPos2f08mfdlpDR64Jj1q3aj2IMjkAjIjraMdijiiEg6S9Jcs/qhWCvP8gWZpidRqyrmMZ7xnURFG6vUisjEcku3Toyep22TN4w91m3mSKzbgu1VzECPdus2yT0oM5Y0y21vQoDPSKqMmVqqgw2VSYaXAISNPggZAnCL+WWYMOwJ7EAexuWDHLFGUCuqMOHXLAXRV5Gz5h1m4pmizGk/IbJucKMXFGGrRjDmHXbWWDWrbO9irE/z8jocdYtEREdJgZ7VDE0yFDhgybJED4JIlt5KxnLt9bsnYxcda4Ce7AH2Ct3jcyeZY+eOes2FkB3NDvr1pHR63RU37ajDh3CXoxhzrotlNEr1l4lAT1246xbIiI7VQIUqe/z+nPfKsRgj4Y4KZvRk6FBhgIfMj4fMjUy/KoEOS3gc9YhOPfzqcgfGBHIvheCnhm0LOMq1lm3gYijvYqzGCMb5KEuV4yhZdurJKLIdGRn3ZbaXsUsxtCKzLotlNHjrFsiIsrHYI+GNH0hUjKzegr8yPj8SEsy/JqMQFqDLy3sRRnOSl0N9oJUo+WKEexFkFvGHQYocbnorFu3pdt2SzFGt1qLRG8UyfYY1Ha/vRijWLDnadatdY8eZ90S0VHG+D/x5bhvFWKwR0OSa5CHANIIIi0HkJKC8AdUSDUKfGEVENnWec6ee24t5qwVutlgT0QBNTvrtjcS0qtu/X0VY+hZPeus255MHIlEVN+j1x6wt1exFmIUnHWrZKtuOeuWiIgGBoM9GpIEZGjwQc0GeWaghxqkEEJAysDn0yDXCPiECp8MSEZGz6jQzUAP9gTsxRnWzF5Yf2nR7KzbWEBfuvWVltHLa6/Spe/RU53tVawZvaKzbhXOuiUi6gsze54w2KMhp1BGL4UapFCDgJSGHxn4ZRX+oAq/nIGQAL8sIGczdpK1MMMa7MmAsFTqihAgwhKUSHbWbcQ567ZARk/U2dqrdKu1SCSMWbch+x69Yku33QC6S5l165bR46xbIjpKGR0YynHfKsRgj4aUXJDns2TzcoFeEiE90IMCv6RC9mmAJBCAiqBfgc+vwRcQ8AVhb8eiIddIWQZEABB+QKmRoYT8SNa4zbotUoxhba+SjCPRG0W6I5w/69Yt0LPOuu1FibNu3ZZuWYxBRER9Y7BHQ4ae0ZOgQrbv0csL9hQEoMAnqZB8AvABqi8FERQIBgA5oAFBS9GGsWfP0l9P+AEtAChBHxLBIHoDYUsfPXt7FevM23bH0m2XqEUyGUOyIwrlUADwUoxhzro1MnpeZ91W6XoDEVFfjJWbcty3CjHYoyFBCAmakPSsnqSHc8ZePWugF8hm9XxQASkbIAoJiuyDKvmhBDJQoMLv1yBpIpv8EpCyCTAhSxCyBMUnQ/XJSPqDSPjD6JEj6EIcnYi7zLjVgzxbMYaazeil9Fm32iEf4FaM4Rbome1VVEDJeJh1a7RX4axbIiIqHYM9GnRCSNmXDE0y9un5LVk9/RVEDQLI6IEerOPTZKjCD0UKIONPIyMr8AcU/V2RfUHogaEkQ5NkZCQ/FMmPpBRCrxxGD6K2YK/PYgwljkRPbtatVmqgZ5t1axRj9CB/6datGIOTMYiIAJRvcaNKt0Az2KPBJfRgT1V9UFQ/MnIAGV8QGcvybRAhJJGBHyp8UM3ATb9cL+bISAHUIIWMFEBa1rN/Ui4chJQNjvQ75PYDJhFCLyJmsNdly+y5F2PYZ90GgfYSl267YZl1q1pm3Xptr8JAj4iISsdgjwZPNm7RFBlKxo9MOoB0IIi0Tw/Cgkhblm0VyNDMoM1gFHQEkUYGgez5emCoh4G5QE9AcgR7ASQRRgJ6Zk8vzIgXzegZs25TnRFopRRjWAM9I4ZTvcy6ZXsVIqI8bL3iCYM9OvKE5aUBWkaC0isj7fcjLQWQCtYgiLTeZsUI9oQKqAKaJiBkkW2hIpmBWw1SyCBgVur6oELOduszGIGhPdgLIZHN7OWKM/ox67bPYgyRK8bgrFsiIjqCGOzRkWe0QlEApADxziGIj/dAjfQi/fkwEqeHEcgu5BqBm9qtYv9fMkj8JYWRn6nBcacB9Y0+pBFADdLmkq+1gMPI7Bn0Sl+/WembQQAJW7BXi07Uep91W0pGL6UBCmfdEhENCPbZ84TBHh1ZloweFOixzzuHoP3571BDPUjH/x3J04/J7thLZzN0KhLdCvZuSWHP2gRO+HoNQiN9CDYGsuFd2lzyNQo49GVffUefQc/s6cFeGsG8YM9Yru171m0/ijHSWpFZt24ZPc66JSKigcFgj44AAfR2AAc6AZ8PiNYBw2J6oJcC4KsB6uqghmuQ8ceRTIUQ8GUQ8GUgS5qeofMLKPUK5ONVqCMySIQUdCOt5//UFBIfZpD4SEE4msHwMQrix+Yye0ZeTMtO5NCDPT13aC3Q6DAbJVtm3abjSCSjSHeGCs+6LViMgRJn3VoDvTQ465aIqA/cs+cJgz0qPwGg42Pgk11AJgLUngA0xHKxzch64MsnQcQzyIwOQyRCCNRk4JcUvWkygGBMQmwyEB7px7CRKaj/lkI3UsggDZ8SwIc7Vez+g4Lhn1Lw2XMz+NSxuWpc/Q6SOX4tF+zpBRq9CKMHMcsMXMus2+4SZt0Wyuhx1i0REQ0BDPaovIyYJdkDdB8AwjGgJ5nL6qUADKsFxtZC1GlQ6xRovRkkkYHPp0KShN48OSwhfJKM8En+bCOWJLqzZRw+LYh9H6v4x5squtIKhndkUAtrsCeZlbhuwV4CYXSLqGPWbRyJRKx/s257jFm3Wj9n3TLQIyIqipk9TxjsUfkYe/NUAJF/A2KnAvVBIHBMLr5JAKjJvoISREAG/AFk5Bok/XpPPXPEmSNoyyCAIILw+9OINKsYqymINajAWAUdsC7jGmPYjGAv4FqNa511m+yN9T3rtgPuS7fmrNti7VWKzbploEdERAOHwR6VjzXYi/0bUNcAHCMBfjm3alkDPe4JZl8BH4TfB8Vfg0RAQEgwZ9pq2U57ete9gNlXL+DPIHy6gjEnKwjIKqSgmg32RF6wV7hAI5I/6/ZgoPSGycas2wSy7VXcRqBx1i0R0YBgZs8TBntUBiL3VQOgSoDmA4Qv9wtqDfbMQE/KZvEENJ8fkGuQ1gBJk6BpMhSfD4rPh4ykL98GEdSrcKUM/EEV/qCCDFT4snW3epyYC/bcCjQSIoSEiKBXjaAnFbPPum33UIxhzLpVjfYqnHVLRFQ2DPY8YbBH5WPN7CnQVyozyMU5RsAXyL78MJdsNdkHIUmAkKAJGYrwQanxQ/EFkEEKQQTM9iy5dit6oGf02JPMxyi0Zy8b7GlhJNIRpLqjSHUd7qzbBArv0WN7FSIiOvIY7FGZZKc+aBKQkXKBnpH4sgZ6ftgCPcgSIOmLsKoAtGxmT9N8UIUfihxARg4gLelZvYCUMQM9WdK/WseqmcGe8OlLwMKPtAgiowWQUkNIZsJIJ0PIdIagHPI66xZ6Ru+wZt0SEZEnRgKhHPetQgz2qAyy6TyhAqovF+QZcU7Q8nIGezL0PXrZrKBQZUDxQ1UkICNDzfihBAJI+wMI+NPw+xQEfNnxaJIKObuzzzk5Q0CGChmK5oeiBZDJBJBWAlDSNVBSQSg9wf5NxuCsWyIiGuIY7NEAsg69NYK9bHWFNeCryX51C/SMYM8YB6vKEKoMofqgZfxARkOmJgN/TQC+YBD+gF5yIUsqfLIGWWhmy5XcU2X37Gk+KKofiupHJhVAJh2ElvBD9AaALg+BHmfdEhENLhXl2V/HPXtEhZiRGXIb9LLBnubTT7EGe0nkgjzr8q0zs5e350/P7omUH2pQggj4oQU1KH4NsqxBkjRIktC/WoM9IUMYS8GKD5oiQ037oKV9EL0y0CvpAVw77AEeZ90SEVEVYLBHh8ktKjNeKqBp+tvWYM9RjAEfCi7j2m6ZkYC0BBGUoAb80AICCABSQACy3nxZkjVAEvrLoEl6hlCRgQwgsnsIRUbSA70e6PFaB/IDPs66JSIaeliN6wmDPRoA1uVbI8tn7EvTACH0P6YB+CX78q1bkJddBTZvoQBQLEUeQf0ewo/c12ywp391BHuqBCiy/R5G8NmLXJxmBHZ9zrpVAM1tMobbrFsjo8diDCIiGhwM9mgAWdNxRpSWDXCMYM+HXKuVJOx79YD8WFF13C5juT6Qvd4vcsGiJOlfzcYr0IM94/o07FXB1k4pRnBXbNZtpj+zbo2/BwZ7REQDgpk9Txjs0QAptJybze6p2VEYPuhxkB/2ggy3fXq2JVzkArS8zKBkv4810BPIDxatvf4SyGX3OlBgBJoAkppejMFZt0REVGEY7NEAcIvQrGul2bVWzacvpRoZPmuwZ71NsWDPuefPuQxsifPyloTdgr0kctk9t0CvF9lZt9b2KsaLs26JiAYFM3ueMNijw+Rot5IX6AX1Pwu9SAIZ5AIz6/JtKbcy2rY4Az23YK9QljCN/MpgI4Yb0Fm3LMYgIqKhgcEe9ZPbsm1eVQVy664y9Dlo2YkakuTeZsX4ai3qNQszkMsKWjN6pQR7bpm9NHKZvV7oQV4ncquynHVLRDQ0Gf9OL8d9q5Dc9ylEhRRqueKMrIyXCmgCUIQ9q2asghp753odLyO+6kIu69ZpeTn32h3Oy2iZlwCQVgA1AXs5rjWjZ126NTJ67KVHRFTN0uk0li5ditGjRyMSiWDy5Ml44YUXSrq2o6MDl156KUaMGIFYLIZp06bhtddeyztPCIF7770Xp59+OuLxOEaOHImZM2diy5Yt/Xpmz8FeT08PbrjhBpxzzjmor6+HLMtYs2ZN3nmXXHIJZFnOe5100kmu933ggQdw0kknIRwOo6mpCatWrfL+09AR5Gy34szqZWDfbJfJHTMCvozQ98NZAz4j6DOCPOv2OOPVBfeA73ACP2s8l1SBdBpQk4DodXxz66zbBPJn3XLploio7ArlGA73VcKevfnz5+OOO+7A3Llzceedd8Lv92PmzJnYvHlz0euEEJg5cyYee+wxLF68GCtWrMC//vUvTJ06Fe+9957t3B/+8Ie4/PLLcdppp+EXv/gFfvjDH2LXrl346le/iu3bt5f2d2TheRn3wIEDuOmmmzB27Fg0Nzdjw4YNBc8NhUJ44IEHIEQuy1FXV5d33n333YfLLrsMc+bMwdVXX40///nPWLx4MRKJBK655hqvj0hHlLOqwgjwfLCvuVo36EmAkHM1DG4vR31HXrsVGfYCD2f7Fq/LuNaAk7NuiYiGtkEq0Ni2bRsef/xx3H777bjqqqsAAHPnzsUpp5yCJUuWYNOmTQWvfeKJJ7BlyxY8+eSTOP/88wEAc+bMQVNTE2644QasXbtWfwRVxb333ouWlha0traa18+ePRvHH388HnnkEXz+85/39GN5DvZGjRqFffv2YcSIEXj11VdxxhlnFL6534+LLrqo6P2SySSuv/56fOMb38Djjz8OAFiwYAFUVcVNN92ESy+91DVApKHCda4ZckGeH3pUZY3MssdV6C1ZhMjeRsoFewHk99SzVuFKyA/2jP16xj2cBRrWFWVrsJcUlm4pXmbdGi/OuiUiOhqsW7cOfr8fCxcuNI/V1NRgwYIF+PGPf4w9e/Zg9OjRrtc++eSTGDlypBnoAUBDQwNaWlrwyCOPIJPJIBAIIJPJIJFIYMSIEbbrjz32WMiyjEgk4vm5PS/jBgKBvAcoRtM0dHV1FXz/pZdewsGDB3H55Zfbjl9xxRXo7u7G+vXrvT4iHRHWjJ4zs1eoBYvxcsyKVYV+irVYwm3fnrPjSaGlXeefOx3Hre/3iNz4MzUJeyluX7Nure1ViIjoiHHbFj4Qrz6yhTt37kRTUxNisZjt+KRJk8z3C3nttdcwYcKEvOOTJk1Cb28vdu3aBUBfFf3CF76A1tZWPProo/joo4/w17/+FRdffDHq6+ttgWapylqg0dvbi9raWtTV1aG+vh6LFi1CT0+P7RxjY+LEiRNtxydOnAhZll03LtJQ4BbsOctonXv2nMGepWgjYynacAv2rIUazpe1gKNQwOcM8ow/d2e/pxnsWd+wZvSM6RgsxiAiOlq1tbWhsbEx73hjYyOEENi7d2+/rgVgu/aRRx5BU1MTvvOd75jb5nbu3IlNmzbh05/+tOfnLlvrlVGjRmHJkiWYMGECNE3D888/j3vuuQd//etfsWHDBsiyHme2tbXB5/OhoaHBdn0gEEB9fX3RvzgabM55uG5799w210mw91rxA8KXLdyQAE3OxY3GUq61p57bEq71kQot45oxaDa4VASgqoBQoQdxbrNunRm9tOVGDPKIiAaF8e/1cty3iEQigZqamrzjoVDIfL8/1wohbNfGYjGcfPLJOPPMM3HWWWdh3759WL58Ob75zW9i06ZNGD58eIk/kK5swd4tt9xi++eWlhaccMIJuP7667Fu3Tq0tLQA0H/4YDDoeo9QKFT0L44GU6FAzxnwycgP+KzXGx2ULVM2jNvKyAV8fpfbOG9p3NYa6FmDPeOlAdA0QFMBYWQajTEanHVLRHRU+Muvgb/+2n4s2VH0knA4jFQqlXc8mUya7/fnWkmSzGtVVcXZZ5+N//iP/8B///d/m+edddZZOPnkk7FixQr87Gc/K/qcTke0qfJVV12FZcuW4YUXXjCDvXA4jHQ67Xp+Mpks+hdHg83ZfsUt0JMcL7icH8jdxyjWAPTiDQX67FsFlsSgYxauW2ZPiFzQp0LfF2gcFwK5DRrWjYLWjB5n3RIRDVkDUY178kX6y2rvDuC+ie7nQ19ydVtxbGtrA6Cvaha71jiv2LUbN27E66+/jl/84he288aNG4cTTzwR//d//1fwexQ7kofaAAAeSUlEQVRyRIO9UCiE+vp6HDx40DzW2NgIVVVx4MAB21JuJpPBJ598UvQvTvc8gJDj2CkATh2ox6ainNW4RjrOmnozSmEB99SbNUg0Kniz1wpZX9YVsj1edMaPbo9kfjsNEJq+XCuM9J5zXloCpRdjMNAjoqPF3wC87jiWHIwHGRKMlnPd3d22Io2tW7dCkiQ0NzcXvdatNcvWrVsRiUTQ1NQEANi/fz8kSYKq5q8pZzIZKIr3KPeITtDo7u7GgQMHcOyxx5rHmpubIYTIaxL4yiuvQNO0on9xuhkALnK8GOgdWW4N7ZwdL41+J9YAy/ln6z8bRR3ZAo5Cs3KTLi8jCWfewliyNZ7BaJRsLcRwFmMkLN/AUkzCQI+IjiqnIv+/sTMG9YkADFpT5dmzZ0NRFNx///3msXQ6jdbWVkyePNlsu7Jv3z688847toBt9uzZ2L9/P5566inz2IEDB7Bu3TrMmjULgUAAANDU1AQhBB577DHb996xYwfeeecd14revpQls5dKpZDJZPJKk2+88UYAwDnnnGMemzZtGoYPH47Vq1djxozc/4BWr16NaDSKc889txyPSIfNCPAkx1dr0Oc2rNZtg50zs2ddBjYyfT4UTuW5PZv1ZfwGWyuDrcFlwvGVs26JiCjfpEmTMGfOHFx33XXYv38/xo0bh9bWVuzevRsPPviged61116LNWvW4IMPPsCYMWMA6MHeHXfcgUsuuQRvvPEGGhoacM8990DTNPz0pz81r50wYQK+9rWv4aGHHkJHRwe+/vWvY+/evVi1ahWi0SiuvPJKz8/dr2Dv7rvvRnt7O/bs2QMAeOaZZ/DRRx8BABYvXoyDBw/i9NNPx0UXXYTPfe5zAIDnn38ezz33HGbOnIlZs2aZ9wqFQrjpppuwaNEitLS0YPr06di4cSMeffRR3HrrrRg2bFh/HpHKzhq8qbAHe5rlmHGuwblvz6iYMP5Zgb0awxns9ZWMti4rW+9pvNKOl3V0hjPjaGTz2DCZiGhIGaQJGgDw8MMPY9myZVi7di0OHTqE8ePHY/369ZgyZYp5jiRJZtcRgyzLeO6553DNNdfgrrvuQiKRwKRJk7BmzRqccMIJtnOfeeYZrFy5Eo899hj++Mc/IhgM4itf+QpuvPHGvHNLIQnrLLMS/fu//zs+/PBD1/fef/991NXVYfHixdi6dSv27t0LVVUxbtw4fOc738HVV18Nn8+Xd90DDzyA22+/He+//z6OO+44fO9738P3vve9gs+wY8eObG++SwHk962hI8WYiOGDvT+K3/Fyzjsr9p4xD83nuG+pwV6x3n/Wnn9GgYaz6tbaYoV79IiI7NoA3I9XX321X0uKh8P8b/8lrwIjy/C99+0AHpw4KD9bOfUrs/f+++/3ec5DDz3k6Z4LFizAggUL+vM4NOisWbq+AqNsXz1bMOYs8PC5vKx9Voot5Tozes4me0aLdGuDZ866JSKqKMa/zstx3yp0RKtxqVpZ++W5LdlazzMCPD/s++qcffmsGb1CHZQLPYszo+e2k9e6d8+tvQqXbomIqDow2KPDZN0fZwRiMnIBV7FrrC8j22cEdD7kMn2lZvWs93bO6bWW8zr37hlBHituiYgqwiBN0KhUDPZoADiXO43xF5rlfeu5haZuWLN4boEe0Hew57ZXzzpGwzmn17p/r1iASkREQ8YgFmhUIgZ7VAbWtixGxk+1vGfN/FmDP7e9eW4z0Ur53s5gzzps19qsz8j0MaNHRETVicEelYG1/Yls+bPxnuw47oN7sAdYZqR5+N7W4M46ENetM7PxPosxiIgqBjN7njDYowFmDfSs/2wN/oyXdbnXKNJwC/ZKbaZsfL9CwZ412+fM6DHQIyKi6sRgjwZQoaDJGuS5BX/WDJ9zGdfr1Ay3ebvOyR7GV866JSKqSGy94gmDPRpARtDkVuTgzOpZs3vGe86ijFKXb51BXrEiEOc+PiIiourGYI/KxFiWNbj1xzPGrbll+7wu3zoDPbe5u87xbEREVJHYesUTBntUBqUsjRqBHpCfyTMCw1IDPcCeqXMu6Vr37XHZloiIji4M9ugIcRZuGHv1rO8ZAWB/llfdAjlnRo9BHhFRVWA1ricM9ugIsfbec5Is53hZurXe2xpIWr+fdXmXiIjo6MNgj44gaybPGpiVOh2jmEKZPRZhEBFVHWb2PGGwR0eYW4bPKNIAvDdQdr6c7xMRER3dGOzRIChUJGEdr1bqfZi9IyI66rDPnicM9mgIcdt7V8r5REREVAiDPRpiChVxuGGgR0R0VDKaLJTjvlWIwR4NQQziiIioCGMYUjnuW4VK3RxFRERERBWImT0iIiKqLGy94gkze0RERERVjJk9IiIiqixsveIJM3tEREREVYyZPSIiIqosbL3iCTN7RERERFWMmT0iIiKqLKzG9YSZPSIiIqIqxsweERERVRZO0PCEmT0iIiKiKsbMHhEREVUW9tnzhJk9IiIioirGzB4RERFVFvbZ84TBHhEREVUWtl7xhMu4RERERFWMmT0iIiKqLGy94gkze0RERERVjJk9IiIiqixsveIJM3tEREREVYyZPSIiIqosbL3iCTN7RERERFWMmT0iIiKqLOyz5wkze0RERERVjJk9IiIiqizM7HnCzB4RERFRFWNmj4iIiCpLufrhsc8eEREREVUaBntERERUWdQyvvqQTqexdOlSjB49GpFIBJMnT8YLL7xQ0mN3dHTg0ksvxYgRIxCLxTBt2jS89tprfV4zYsQIyLKMp556qqTv48Rgj4iIiCqLUaAx0K8Sgr358+fjjjvuwNy5c3HnnXfC7/dj5syZ2Lx5c9HrhBCYOXMmHnvsMSxevBgrVqzAv/71L0ydOhXvvfdeweuWLVuGZDIJSZL6frgCGOwRERERlWDbtm14/PHHsXz5cixfvhzf/e538eKLL2Ls2LFYsmRJ0WufeOIJbNmyBQ899BCuv/56XHbZZXjppZfg8/lwww03uF7z+uuv495778XSpUsP67kZ7BEREVFlGaTM3rp16+D3+7Fw4ULzWE1NDRYsWIAtW7Zgz549Ba998sknMXLkSJx//vnmsYaGBrS0tOB3v/sdMplM3jVXXnklLrjgAnzpS1+CEKL4wxXBYI+IiIioBDt37kRTUxNisZjt+KRJk8z3C3nttdcwYcKEvOOTJk1Cb28vdu3aZTv+xBNPYOvWrfj5z39+2M/NYI+IiIgqiwIgU4ZXH61X2tra0NjYmHe8sbERQgjs3bu3X9cCsF2bTCZxzTXX4Ac/+AGOO+644g9VAgZ7RERERCVIJBKoqanJOx4Khcz3+3OtEMJ27c9+9jMoioLrrrtuAJ6aTZWJiIio0qgA+l+cqhO/BvBrx8GOopeEw2GkUqm848lk0ny/P9dKkmRe+8EHH2DlypVYvXo1IpFI8Z+hRAz2iIiI6OgjXQTgIvsxsQPAxIKXNDY2ui7VtrW1AQBGjRpV9FrjvGLX/uQnP8GnPvUpfOUrX8Hu3btt5/zrX//C7t27MWbMGE+tWBjsERERUeXpf3FqvzU3N2PDhg3o7u62FWls3boVkiShubm56LWbNm3KO75161ZEIhE0NTUBAD766CO8++67OP74423nSZKEyy67DJIk4dChQ6itrS35ublnj4iIiKgEs2fPhqIouP/++81j6XQara2tmDx5MkaPHg0A2LdvH9555x2oqmq7dv/+/bYpGAcOHMC6deswa9YsBAIBAMAtt9yCp59+Gr/97W/N18033wwAWLp0KZ5++mlEo1FPz83MHhEREVEJJk2ahDlz5uC6667D/v37MW7cOLS2tmL37t148MEHzfOuvfZarFmzBh988AHGjBkDQA/27rjjDlxyySV444030NDQgHvuuQeapuGnP/2pee2ZZ56Z933r6uoghMAZZ5yBWbNmeX5uBntEREREJXr44YexbNkyrF27FocOHcL48eOxfv16TJkyxTxHkiTIsn3xVJZlPPfcc7jmmmtw1113IZFIYNKkSVizZg1OOOGEPr/v4YxLk8ThtGQeRDt27MDEiRMBXAogv28NERERlUMbgPvx6quvujYJLqfcf/tfBVCO760XaAzGz1ZO3LNHREREVMUY7BERERFVMe7ZIyIiogpjzEsrx32rDzN7RERERFWMmT0iIiKqMArKk4VjZo+IiIiIKgwze0RERFRhuGfPC2b2iIiIiKoYM3tERERUYVSUJwun9n1KBWJmj4iIiKiKMbNHREREFYZ79rxgZo+IiIioijGzR0RERBWGmT0vmNkjIiIiqmLM7BEREVGFYTWuF8zsEREREVUxZvaIiIiownDPnhcM9oiIiKjCcBnXCy7jEhEREVUxZvaIiIiownAZ1wtm9oiIiIiqGDN7REREVGEUlCcLx8weEREREVUYZvaIiIiownDPnhfM7BERERFVMWb2iIiIqMKwz54XzOwRERERVTFm9oiIiKjCcM+eF8zsEREREVUxZvaIiIiownDPnhfM7BERERFVMWb2iIiIqMJwz54XDPaIiIiownAZ1wsu4xIRERFVMWb2iIiIqMJwGdcLZvaIiIiIqhgze0RERFRhmNnzgpk9IiIioirGzB4RERFVGAXlycIxs0dEREREFYaZPSIiIqow3LPnBTN7RERERFWMmT0iIiKqMJyg4QUze0RERERVjJk9IiIiqjDcs+cFM3tEREREVcxTsLd9+3YsWrQIp5xyCmKxGMaOHYsLL7wQf//73/POffvttzFjxgzE43HU19dj3rx5OHDggOt9H3jgAZx00kkIh8NoamrCqlWr+vfTEBER0VHA2LM30C/u2cNtt92Gp59+GmeffTbuvPNO/L//9/+wceNGTJgwAW+++aZ53p49e/DlL38Z//jHP7B8+XJcc801WL9+Pb7+9a9DUewp0vvuuw8LFy7EqaeeilWrVuHMM8/E4sWLsWLFioH5CYmIiKjKGMu4A/3qexk3nU5j6dKlGD16NCKRCCZPnowXXnihpKfu6OjApZdeihEjRiAWi2HatGl47bXXXM/dvHkzvvSlLyEajaKxsRFXXnklenp6Svo+Tp727F199dX49a9/Db8/d1lLSwtOPfVULF++HGvWrAEA3HLLLUgkEti5cydGjx4NADjjjDPwta99Da2trfjud78LAEgmk7j++uvxjW98A48//jgAYMGCBVBVFTfddBMuvfRS1NXVlfBkfwNwqpcfhSoGP9vqxc+2uvHzrV7vDvYDDKr58+fjqaeewlVXXYVx48ahtbUVM2fOxIYNG3DmmWcWvE4IgZkzZ+Jvf/sblixZgvr6etxzzz2YOnUqduzYgc985jPmuTt37sTZZ5+Nk046Cb/4xS/wz3/+EytWrMC7776L9evXe35mT5m9yZMn2wI9ABg3bhxOPvlkvPXWW+axp556Cuedd54Z6AHAWWedhaamJvzmN78xj7300ks4ePAgLr/8cts9r7jiCnR3d3v4gV738mNQReFnW7342VY3fr7VaygEe4OzjLtt2zY8/vjjWL58OZYvX47vfve7ePHFFzF27FgsWbKk6LVPPPEEtmzZgoceegjXX389LrvsMrz00kvw+Xy44YYbbOf+6Ec/wvDhw/GnP/0Jl156KW688UasWrUKzz//fMlZRKsBKdDYv38/GhoaAAB79+7Fxx9/jM9//vN5502aNMmWrjT+PHHiRNt5EydOhCzLBVObREREREfaunXr4Pf7sXDhQvNYTU0NFixYgC1btmDPnj0Fr33yyScxcuRInH/++eaxhoYGtLS04He/+x0yGb26uKurCy+88ALmzp2LaDRqnjtv3jxEo1Fb0qxUhx3srV27Fnv27MF//ud/AgDa2toAAI2NjXnnNjY24uDBg+YP1NbWBp/PZwaKhkAggPr6euzdu/dwH4+IiIiqzuDs2du5cyeampoQi8VsxydNmmS+X8hrr72GCRMm5B2fNGkSent7sWvXLgDA3/72NyiKkpcICwQCaG5u7lci7LCCvbfffhuLFi3ClClTMG/ePABAIpEAoEe6TqFQyHZOIpFAMBh0vXcoFDLPIyIiIhpsbW1tBZNZQoiiSapi1wIwr21ra4MkSQXP7U8irN9Nlffv349zzz0XxxxzDJ544glIkgQACIfDAIBUKpV3TTKZtJ0TDoeRTqdd759MJs3z3OQCwWcAdAB40PLuuOyLKl8SQNtgPwSVBT/b6sbPtzq8i/w9evsBYJATMvtQngbI7i3iDIlEoqRkltdrhRC2RBhQOGnWn7/3fgV7nZ2dmDFjBjo7O7Fp0yaMHDnSfM+IRI3lXKu2tjYMHz4cgUDAPFdVVRw4cMC2lJvJZPDJJ59g1KhRBZ/hgw8+yP5pX/brh5Z3PwTwv/34yWhoun+wH4DKhp9tdePnW80++OADTJky5Yh+z4aGBkQiEfT2PlW27xEMBvO2lxnC4XBJySyv10qSZEuEAYWTZsW+RyGeg71UKoXzzjsP7777Ll588UV89rOftb0/atQoHHvssdi+fXvetdu2bUNzc7P5z83NzRBCYPv27ZgxY4Z5/JVXXoGmabZznaZPn461a9fi05/+dL9+cCIiIvIukUjggw8+wPTp04/49x4zZgzeeuutgkMaBkJDQwPGjBnj+l6hZVQjwVUsSdXY2FgwEWa91lgSLnRuse9RiKdgT9M0tLS04OWXX8Yzzzxjbkh0uuCCC7BmzRrs2bPHbL/y4osvYteuXbj66qvN86ZNm4bhw4dj9erVtmBv9erViEajOPfccws+S0NDA7797W97eXwiIiIaAEc6o2c1ZsyYgsFYuTU3N2PDhg3o7u62FWls3boVkiQVTVI1Nzdj06ZNece3bt2KSCSCpqYmAMApp5wCv9+P7du3Y/bs2eZ5mUwGO3fuxIUXXuj5uT0VaPzgBz/A73//e5xzzjk4cOAAHnnkEdvL8KMf/QiRSARTp07FqlWr8LOf/QwtLS047bTTcPHFF5vnhUIh3HTTTXj22WfR0tKCBx54APPnz8ejjz6K66+/HsOGDfP8AxERERGVw+zZs6EoCu6/P7dFIZ1Oo7W1FZMnTzYTXPv27cM777wDVVVt1+7fvx9PPZVbgj5w4ADWrVuHWbNmmVvcamtrcfbZZ2Pt2rW2iRlr1qxBT08PWlpavD+48GDq1KlCluWCL6s333xTzJgxQ8RiMTF8+HAxb9488fHHH7ve91e/+pU48cQTRSgUEieccIK48847vTwWERER0RHR0tIigsGgWLJkibj//vvFmWeeKYLBoNi0aZN5zvz584UkSWL37t3mMVVVxRe/+EVRW1srbrzxRnHPPfeIU045RdTV1Yldu3bZvseOHTtEOBwWEyZMEPfee6/48Y9/LMLhsDjnnHP69cySEEJ4DxGJiIiIjj7pdBrLli3D2rVrcejQIYwfPx4333wzzj77bPOcSy65BA8//DD+8Y9/2JacOzo6cM011+C3v/0tEokEJk2ahJUrV+L000/P+z6bN2/G0qVLsWPHDsTjcVx44YW49dZbbY2WS8Vgj4iIiKiKDci4tCMtnU5j6dKlGD16NCKRCCZPntyvWXE0uP70pz9BluW8l8/nw7Zt22znvv3225gxYwbi8Tjq6+sxb968slZjUel6enpwww034JxzzkF9fT1kWcaaNWtcz/XyOT7wwAM46aSTEA6H0dTUhFWrVpXzx6ACSv18L7nkEtff55NOOsn1vvx8B9/27duxaNEinHLKKYjFYhg7diwuvPBC/P3vf887l7+7la3fTZUH0/z58/HUU0/hqquuwrhx49Da2oqZM2diw4YNOPPMMwf78cij73//+3mzlMeNyzXF3rNnD7785S/jmGOOwfLly9HV1YUVK1bg9ddfx7Zt2+D3V+T/jKvGgQMHcNNNN2Hs2LFmpZobL5/jfffdh8suuwxz5szB1VdfjT//+c9YvHgxEokErrnmmiP0kxFQ+ucL6EV3DzzwAKwLRnV1dXnn8fMdGm677TZs3rwZc+bMwfjx47Fv3z7cddddmDBhAl5++WUzUOfvbhXo106/QfTyyy8LSZLEf/3Xf5nHksmkGDdunJgyZcogPhl5tWHDBiFJknjyySeLnnfZZZeJaDQq/vnPf5rHXnjhBSFJkvjlL39Z7sekPqTTabF//34hhBDbt28XkiSJhx56KO+8Uj/HRCIhGhoaxKxZs2zXf+c73xHxeFy0t7eX6SchN6V+vhdffLGIx+N93o+f79CxZcsWkclkbMf+/ve/i1AoJObOnWse4+9u5au4Zdx169bB7/dj4cKF5rGamhosWLAAW7ZswZ49ewbx6ai/uru7bSXqVk899RTOO+88s6QdAM466yw0NTXhN7/5zZF6RCogEAhgxIgRfZ5X6uf40ksv4eDBg7j88stt119xxRXo7u7G+vXrB+7hqU+lfr4GTdPQ1dVV8H1+vkPH5MmT81ZGxo0bh5NPPhlvvfWWeYy/u5Wv4oK9nTt3oqmpydbMEIDZ4Hnnzp2D8Vh0GC655BLU1tYiFAph2rRpePXVV8339u7di48//jhvmRfQP/PXXnvtSD4q9ZOXz9H488SJE23nTZw4EbIs8zMfwnp7e1FbW4u6ujrU19dj0aJFtj5hAD/fSrB//35zXBh/d6tDxW12amtrM+fvWhnjRdzGmNDQFAwGMXv2bMycORMNDQ148803sXLlSnzlK1/B5s2bcdppp5njYgp95gcPHkQmkzGbUdLQ5OVzbGtrg8/ny5tNGQgEUF9fz9/xIWrUqFFYsmQJJkyYAE3T8Pzzz+Oee+7BX//6V2zYsAGyrOcW+PkObWvXrsWePXtw8803A+DvbrWouGAvkUigpqYm73goFDLfp8rwxS9+EV/84hfNfz7vvPNwwQUXYPz48bjuuuvwhz/8wfw8+/rMGewNbV4+x0QigWAw6HqfUCjE3/Eh6pZbbrH9c0tLC0444QRcf/31WLdundn1n5/v0PX2229j0aJFmDJlCubNmweAv7vVouKWccPhMFKpVN7xZDJpvk+V6zOf+Qy++c1v4qWXXoIQwvw8+ZlXNi+fYzgcRjqddr1PMpnk511BrrrqKkiSZGuNxc93aNq/fz/OPfdcHHPMMXjiiScgSRIA/u5Wi4oL9hobG820spVxbNSoUUf6kWiAHXfccUin0+jp6TGXDgp95sOHD2dWrwJ4+RwbGxuhqmpeD69MJoNPPvmEv+MVJBQKob6+HgcPHjSP8fMdejo7OzFjxgx0dnbi+eefx8iRI833+LtbHSou2GtubsauXbvQ3d1tO75161ZIkoTm5uZBejIaKO+99x5CoRBisRhGjRqFY489Ftu3b887b9u2bfy8K4SXz7G5uRlCiLxzX3nlFWiaxs+8gnR3d+PAgQM49thjzWP8fIeWVCqF8847D++++y7Wr1+Pz372s7b3+btbHSou2Js9ezYURcH9999vHkun02htbcXkyZNtpeE0tLl1X//LX/6C3//+95g+fbp57IILLsCzzz5ra6vz4osvYteuXeY+IBr6Sv0cp02bhuHDh2P16tW261evXo1oNIpzzz33iD0zlSaVSuX9H3AAuPHGGwEA55xzjnmMn+/QoWkaWlpa8PLLL2PdunVmVwsn/u5WvoqcjXvhhRfit7/9Lb7//e+bEzS2b9+O//3f/8WUKVMG+/GoRGeddRbC4TDOPPNMjBgxAm+88QZ++ctfoqamBps3bzb/H+Y///lPTJgwAXV1dbjyyivR1dWFlStXYsyYMdi2bRuXcYeAu+++G+3t7dizZw/uvfdefOtb3zIHey9evBjxeNzT57h69WosWrQIF1xwAaZPn46NGzdi7dq1uPXWW7F06dLB+jGPWn19vgcPHsTpp5+Oiy66CJ/73OcAAM8//zyee+45zJw5E88++6ztfvx8h4bvf//7uPPOOzFr1izMmTMn7/1vf/vbALz9O5if7RA1iA2d+y2VSoklS5aIUaNGiXA4LL7whS+I//mf/xnsxyKP7rrrLjF58mTR0NAggsGgGD16tJg/f75477338s598803xYwZM0QsFhPDhw8X8+bNEx9//PEgPDW5+fSnPy1kWXZ97d692zzPy+f4q1/9Spx44okiFAqJE044Qdx5551H6schh74+3/b2djFv3jzR1NQkYrGYCIfD4tRTTxW33XabUBTF9Z78fAff1KlTC36usizbzuXvbmWryMweEREREZWm4vbsEREREVHpGOwRERERVTEGe0RERERVjMEeERERURVjsEdERERUxRjsEREREVUxBntEREREVYzBHhEREVEVY7BHREREVMUY7BERERFVMQZ7RERERFWMwR4RERFRFWOwR0RERFTF/j8w3i5R8W1TDgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x000000002436B320>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.colorbar.Colorbar object at 0x00000000296C6A58>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imshow(M); colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The contour computation is not 100% robust at this time. In practice this is not a serious issue, since the coordinates of $\\vt{x}$ depend on the location of the e.g. Gauss-Legendre quadrature points used for the testing integration in the variational method this toolbox is meant for. These coordinates are *weird* numbers as in highly unlikely to hit the geometric rountine's boundary cases."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.6",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
